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Abstract 

This article reviews recent developments in multiresolution analysis which make it a pow- 
erful tool for the systematic treatment of the multiple length-scales inherent in the electronic 
structure of matter. Although the article focuses on electronic structure, the advances described 
are useful for non-linear problems in the physical sciences in general. Among the reviewed de- 
velopments is the construction of exact multiresolution representations from extremely limited 
samples of physical fields in real space. This new and profound result is the critical advance in 
finally allowing systematic, all electron calculations to compete in efficiency with state-of-the-art 
electronic structure calculations which depend for their celerity upon freezing the core electronic 
degrees of freedom. This review presents the theory of wavelets from a physical perspective, 
provides a unified and self-contained treatment of non-linear couplings and physical operators 
and introduces a modern framework for effective single-particle theories of quantum mechanics. 

PACS numbers: 71.15.-m, 31.15.-p, 02.30.Mv 
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Introduction 



The focus of this review is the application of wavelet theory to the determination of elec- 
tronic structure. Wavelet theory has at its foundation a single, simple idea: multiresolution 
analysis\ Mal89 , Mey86 , Mey90fl , a relatively recent and mathematically rigorous theory of the 
description of functions which provides simultaneously for a homogeneous underlying description 
of space and the capacity to control and vary the resolution of this description at will. 

Interest in multiresolution analysis and wavelet theory has mushroomed dramatically since their 
introduction. Over seventy-two monographs have been written on the general subject of wavelets 
within the last six years. (For a comprehensive review of the literature of the field prior to 1993, 
see ||PSU93|| .) In addition to several recent introductory texts [ |Kai94| , |HW96| , |SN96j , |GGBB97|] , 
specialized monographs are now available which discuss applications to such wide-ranging fields as 
chemical engineering ||MJ94|1 , bio-medical engineering [|Aka9^1 an d applied science ||bFB97|1, as well 
as the traditional areas of application in mathematics | MCC97|| and signal processing! 3ut98 |. In 
terms of recent monographs, those of the greatest relevance to the present discussion describe the 



application of wavelets methods to partial differential equations [DK097, ALM + 97 |. 

The well-known fact that the electronic wave functions in molecular and condensed-matter 
systems vary much more rapidly near the atomic nuclei than in interatomic regions calls for pre- 
cisely the capabilities of multiresolution analysis. Figure [l] illustrates the multiscale behavior of 
electronic wave functions, using the carbon atom as an example. The curves in the figure show 



the Kohn-Sham orbitals of the atom as computed within the local density approximation [KS65] to 
density functional theory [HK64]. In the immediate vicinity of the nucleus and its strong attractive 
potential, the electrons possess large kinetic energies, as reflected by high spatial frequencies evi- 
dent in the orbitals. In this example, the high-frequency "core" region extends only approximately 
0.5 Bohr radii out from the nucleus in the case of carbon, beyond which the variations in the wave 
functions are quite smooth. Resolving the cusps in the s states of this atom requires a resolution on 



the order of 0.03 Bohr (corresponding to a plane wave cutoff| PTA + 92 | of nearly 160,000 Rydberg). 
To provide this resolution uniformly throughout a computational cell of 8 Bohr on a side would 
require a basis with 16 million coefficients. The vast majority of these basis functions would be 
wasted as they would serve to provide unnecessarily high resolution outside the core region; only 
about sixteen thousand functions would be needed to provide the required resolution uniformly 
throughout the core region defined above. 

Multiresolution analysis allows us to add resolution precisely into the core in a system- 
atic, hierarchical manner. The solid curves in Figure |] come from the earliest reported ap- 
plication of multiresolution analysis to self-consistent electronic structure calculations in three 
dimensions | ACLT95 1 . Despite the high resolution needed in the core, these multiresolution cal- 



culations required fewer that three thousand basis functions and yet produce results nearly indis- 
tinguishable from the output of atomic calculations carried out at essentially infinite resolution by 
exploiting spherical symmetry to produce an effective one-dimensional problem and then using a 
very fine radial grid. 

Figures present similar all-electron calculations for the N2 molecule | Ari95 | . Figure |2| shows 
the lowest energy Kohn-Sham orbital. The la symmetry of the state and the cusps near the nuclei 
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Figure 1: Early electronic structure calculations using multiresolution analysis [ACLT95]: Kohn- 
Sham orbitals of the carbon atom within the local density approximation from standard atomic 
software (diamonds) and multiresolution analysis (curves). 
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Figure 2: la Kohn-Sham orbital of N2 in plane containing both nuclei. 



are clearly visible. Figure || shows results for the bond length and vibrational frequency of the 
molecule. Using the cubic fit in the figure, we find a bond length and a spring constant within 
0.1% and 7% of experiment, respectively. Both these and the preceding atomic calculations were 
carried out in the same 8 Bohr periodic cell at an effective resolution resolution of corresponding 
to 16 million coefficients. Both sets of calculations used the I = 3 interpolating scaling functions of 
the product form discussed in Sec. V-C.4 and the basis restriction strategy described in Sec. QI-A 
with seven levels of resolution. 

The succeeding pages lay out the explicit details of how these and other calculations of elec- 
tronic structure have been performed using multiresolution analysis. Although the emphasis of this 
review is the calculation of electronic structure, the techniques we describe are widely applicable 
to other physical problems involving coupled sets of linear and non-linear partial differential equa- 
tions. Shortly after the initial reports of the above applications to electronic structure calculations, 
independent applications of the ideas of wavelet theory to physical problems described by partial 
differential equations appeared in a variety of areas, including combustion flFS94 , FS97| and fluid 
mechanics |VYP 9 7 , |CD97 |. General model problems have also been explored [BK95, BN96 , BK97 |. 
And, more recently, the solution of Poisson's equation has been studied [GI98, LAE98|. 

The issue of multiple length-scales in electronic structure is not new. It has driven the de- 
velopment of a variety of techniques which are now quite mature, including the linear muffin tin 
orbital (LMTO) method | And75[| , the linearized augmented plane wave (LAPW) method [gm9l , 
the full potential LAPW (FLAPW) method fWKWF8i|1 and the plane wave pseudopotential 
approach I PTA + 92| . The first three methods use one type of basis set inside of a set of spheres 
organized around the nuclei and another type of basis set outside of the spheres. The wave func- 
tions are then matched at the spherical boundaries to determine the solution. The plane wave 
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Figure 3: Energy of N2 as a function of inter-nuclear separation: results of multiresolution analysis 
(crosses), cubic fit and quadratic with experimental curvature and minimum (solid and dashed 
curves, respectively). 
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pseudopotential approach replaces the atomic core with an effective potential manufactured to 
have similar scattering properties. While each of these approaches has had great success, none is 
systematically improvable to complete convergence in a simple, practical manner, and each requires 
great care and expertise in the selection and construction of the atomic spheres or in the devel- 
opment of appropriate pseudopotentials. As a result, a general method is still needed to obtain 
unambiguous results to a sufficient accuracy to permit direct and systematic study of the relative 
accuracy of competing density functionals and alternate theories of electronic structure. 

As illustrated briefly above and discussed in depth below, multiresolution analysis provides 
a systematic approach which can replace millions of grid points with merely thousands of basis 
functions. The development of the multiresolution analysis of electronic structure therefore holds 
the promise of at last enabling the systematic evaluation of different theories of electronic structure 
at high precision. In addition, the mathematical structure of multiresolution analysis is sufficiently 
rich so that new algorithms and techniques are constantly being developed, making it probable 
that multiresolution approaches will prove to be not only more accurate and systematic, but also 
more computationally efficient than present approaches. 

This review is organized as follows. Section [n] overviews the equations of density functional 
theory and gives an extremely brief review of other modern, systematically improvable approaches 



to the calculation of electronic structure. Section III introduces a new basis-set independent ma- 



trix language to express the equations of density functional theory. Section IV lays down the 
mathematical framework of multiresolution analysis in a language suited for physical applications 
in multiple dimensions, and Sec. |V] gives specific examples of basis functions which fit into this 



framework. Sees. VI VII then describe new methods which are needed to make the application 



of the operators and transforms associated with these functions feasible in physical calculations of 



complex systems. Finally, the review concludes in Sec. VIII with a few brief remarks. 

II Electronic Structure 
II-A Kohn-Sham Lagrangian 

Over the last several decades, density functional theory has proven an accurate, reliable and effective 
tool for predicting electronic structure. It has found application in such diverse areas as the study of 
surfaces, point defects, melting, diffusion, plastic deformation, disorder, catalysis, phase transitions 
and chemical reactions. For reviews see |poh84 |AFQ85| , |Pic89| , pTA + 92f| . 



In standard atomic units, fr = m = e = 1, the equations of density functional theory in the local 
density approximation (LDA) |KS65 | are equivalent to finding the saddle point of lowest energy of 



the Lagrangian functional, 

CLDA(mA) = ^E// d3r HV^(r)|| 2 + |d 3 ry ion (rHr) 

+ f d^r e xc (n(r))n(r) d 3 r (ft(r)(n(r) — no) 

|d 3 r||V<A(r)|| 2 , (1) 

where the electron density is defined as 

n(r) = £/hMr)| 2 , (2) 



1 

87 
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the orbitals {ifti} obey the ortho-normality constraint 



d 3 rtft*(r)tft j (r)=6 ij , (3) 

and 3i(z) denotes the real part of the complex number z. Here, Vj on (r) is the potential of each 
electron due to the presence of the nuclei (and core electrons in the case of pseudopotential cal- 
culations) in the system, and e xc (n) is the exchange-correlation energy per electron in a uniform 
electron gas of density n. For simplicity, in this review we hold the occupation factors / in (|l|) 
and (§) fixed at / = 2 to reflect the fact that two electrons of spin 1/2 may occupy each orbital. 
For calculations in periodic systems, we introduce no, which corresponds to a positive charge back- 
ground neutralizing the electronic charge density. The effect of this background on the total energy 
is properly accounted when the Ewald summation is used to compute the interionic interactions. 
At the saddle point, the value of Clda is the total Kohn-Sham energy of the system[ KS65[| , and 



the fields {ifti(r)} and <f>(r) are the Kohn-Sham orbitals (electronic wave functions) and the Hartree 
potential (the electrostatic field arising from the mean electron density), respectively. Taking the 
real part of the integral coupling the electron density n(r) to the Hartree field (j){r) ensures that (ft(r) 
is real at the saddle point. This Lagrangian formalism for density functional theory, introduced 



in [ LAE98 ], has the advantage over the standard energy functional approach of rendering local 
in space all couplings among the physical fields. This not only dramatically simplifies formal 
manipulations but also allows for the practical strategy of a direct search for the saddle point to 
solve the Schrodinger and Poisson equations simultaneously. 

Three factors make locating the saddle point of (|l]) challenging. First, the Lagrangian deals with 
continuous fields which must be describe in terms of a finite number of coefficients for the purpose 
of calculation. Second, the solutions we seek exhibit multiscale behavior. Finally, the Lagrangian 
couples the fields non-linearly, both through the exchange-correlation energy density e xc (n) ■ n and 
through the term coupling the Hartree field and the electronic charge density, (ft(r)n(r). 

A variety of systematically improvable approaches have appeared in the literature to meet these 
challenges. To place the development of multiresolution analysis in context, we now give a brief 
overview of these other approaches. We shall not discuss the muffin-tin families of approaches, 
which are not so closely related to multiresolution analysis, beyond their brief mention in the 
introduction. 

II-B Systematic basis approaches 
II-B.l Plane wave approach 

The plane wave approach is reviewed in detail in [PTA + 92]. In this approach, the Kohn-Sham 



orbitals and the Hartree potential are expanded in a discrete basis of plane waves (complex ex- 
ponentials) consistent with periodic boundary conditions. The resulting discrete set of expansion 
coefficients for the orbitals {ifti} and Hartree potential (ft are well-suited for computation. 

In a plane wave basis, differential operators are diagonal, making their implementation partic- 
ularly simple. The remaining couplings in the Lagrangian (jl]) are the non-linear, spatially local 
couplings. When using a plane wave basis, one implements these couplings on a point by point 
basis in real space, using the fast Fourier transform (FFT) to convert efficiently between the real 
space and the plane wave representations. 

Plane wave calculations may be brought to convergence simply by increasing a single parameter, 
the kinetic energy below which all plane waves are included in the basis. This makes systematic 
basis set convergence studies straightforward. However, the extremely high resolution required 



10 



near atomic nuclei combined with the uniform resolution afforded by plane waves makes the direct 
application of this method prohibitive for all but the lightest elements. The introduction of pseu- 
dopotential theory overcomes this limitation, but at the cost of introducing the pseudopotential 
approximation, which is uncontrolled. One great advantage of multiresolution analysis is that it 
maintains the regularity of plane wave expansions while allowing variable resolutions and thus the 
direct treatment of heavier elements. 

Finally, even the wave functions in pseudopotential calculations at times require significant 
resolution near ionic cores, particularly when dealing with first-row elements or transition metals. 
This opens the exciting possibility of combining the pseudopotential approach with multiresolution 



analysis, an issue which has begun to be explored! WC96j 



II-B.2 Finite element and adaptive mesh approaches 

There is an extensive literature dedicated to the finite element approach, particularly in the fields 



of solid [ Br a97] and fluid [Chu92b| mechanics. See [Whi97| for a review of recent developments, and 
flCB9? 1 for an introduction. 

The application of the finite element approach to electronic structure calculations began in the 
late eighties [ |WWT89 ] and has undergone a recent revival[TT95b, PT95a, TT96]. As in the plane 
wave approach, the finite element method expands the electronic orbitals in terms of a set of basis 
functions. By using localized basis functions concentrated in the regions of space requiring high 
resolution, these bases can provide a much more efficient description of electronic wave functions. 

Finite elements also represent an extremely efficient method for dealing with non-linear interac- 
tions by providing, as do plane waves, highly efficient rapid transforms. For finite elements, these 
transforms are based on two properties of the basis functions, cardinality and interpolation [Whi97|. 
Below, we shall discuss how to construct multiresolution analyses which maintain these two highly 
desirable properties. 

One great difficulty with the application of finite elements to electronic structure is that finite 
elements basis sets must be uprooted and reformed as the nuclei move. Each coefficient in a finite 
element representation corresponds to the value or weight of a function over the region of one 
basis function and thus cannot be taken to be small where the electronic orbitals themselves are 
non-negligible. As a result, coefficients associated with the basis functions which are uprooted 
as the atoms moves carry large values, and this process must be managed with extreme care. 
Multiresolution analysis provides an elegant solution to this difficulty which we discuss briefly 



immediately below and in more depth in Sec. III-A 



Another potential solution to variations in the basis as the atoms move is provided by a branch 
of finite element methods particularly attractive for the calculation of electronic structure, the 
"Ricmannian metric" | Gyg95 , pyg95| , GG95, DKCJ94| or "adaptive curvilinear-coordinate" [ Ham95 , 

TT96 | approach. This approach lays down the finite 



Ham96b , Ham96a . |Ham97 



ZMK9C, 



NC97 



element mesh according to a smooth mapping from a underlying cubic grid of points, thereby 
ameliorating the problems of uprooting the finite element grid as the atoms move. Because it 
preserves the underlying cubic topology of the grid, the Riemannian metric approach falls into 
the class of "structured mesh" methods, for which it is difficult to generate very strongly graded 
meshes [Ran93], a limitation which multiresolution analyses do not suffer. 



II-B.3 Multigrid algorithms 

Multigrid algorithms provide an extremely effective means of solving equations whose convergence 
is limited by a wide range of length-scales. This approach too has an extensive literature associated 
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with its application in a wide variety of fields. For an in-depth introduction, see [Wes92]. 

Explorations of the multigrid approach as applied to electronic structure calculations also 



date back to the late eighties [WWT89] and have become much more common in the last few 
years [ GG95 , BSB96 , TB95| , BBSB97| . Multigrid algorithms do not specify a basis set and leave 
open the issue of how best to discretize physical problems. The mathematical structure of multigrid 
algorithms parallels very closely the ideas of multiresolution analysis, and multigrid algorithms can 
be applied directly or easily generalized to the solution of differential equations expressed in wavelet 
bases HRJZ94 [Yes9"7| , [YA98H . 



II-B.4 Multiresolution analysis 

Wavelet bases place functions of varying resolution on a multiresolution grid while maintaining 
a uniform resolution throughout all of space in the precise mathematical sense of multiresolution 
analysis | Mal8E| , Mcy86 , Mcy9C j . The mathematical regularity of the resulting basis leads to efficient 
fast transforms [Dau9^, Chu92a] and methods to apply differential operators[BCR91, L AE98j ] . 

In contrast to the expansion coefficients of a finite-element expansion which reflect directly the 
values of a function, the coefficients of a multiresolution analysis separate information into different 
length-scales. This subtle but critically important difference means that, so long as a function 
varies smoothly, the fine-scale coefficients will be quite small even where the value of a function 
is quite large. As described in Sec. III-A, this means that, as the atoms move, one may arrange 
for the changes in the basis to involve the truncation of only coefficients which are small, thereby 
effectively providing a high resolution throughout all of space with an extremely limited number of 
coefficients. This also means that no particular care is needed to handle the regriding as the atoms 
move. 

The first electronic structure calculations to use such a basis employed a discrete frame of 



non-orthogonal Gaussian- Mexican hat basis functions [CAJL93|. This work established the effi- 
cacy of such bases for representing electronic wave functions, but the calculations were limited 
to one-electron systems. The first reported self-consistent, multiple-electron density functional 
calculations [ |ACLT95 , Ari95] used the semicardinal bases described in Sec. [V] and employed the an- 
alytically continued conjugate gradient approach [APJ92] to solve the Kohn-Sham equations. How 
Poisson's equation was solved in these calculations is described in more detail in |[LAE98||. Wei and 



Chou[WC96] carried out the first calculations employing orthogonal Daubechies wavelets [Dau92]. 
This work studied molecular dimers within the local density approximation using Daubechies D6 
wavelets, employed self-consistent iteration to solve of the Kohn-Sham equations, and represented 
the first use of pseudopotentials in the multiresolution analysis of electronic structure. Since that 



time, Tymczak and Wang|TW97] used Daubechies D8 wavelets and introduced the innovations of 
dynamically refining the basis and the use of the Car-Parrinello approach to solve the Kohn-Sham 
equations. Most recently, Goedecker and Ivanov | GI9Sfl have applied lifted wavelets [3we96] to the 
solution of Poisson's equation. 

Until recently, the primary bottleneck in multiresolution analysis calculations of electronic struc- 
ture had involved the performance of transforms and the application of differential operators. The 
standard transforms associated with orthogonal wavelet bases require or produce the values of 
the electronic orbitals on a uniform grid at the finest resolution. To use these transforms, one 
"unpacks" each electronic orbital from its stored coefficients, operates on the unpacked version, 
and then "repacks" the result. Although this gives great benefit in terms of the use of memory, 
processing the wave functions in their highly redundant unpacked representation still involves sig- 
nificant memory and also much wasted computation. For example, applying this approach to the 



12 



nitrogen dimer calculations described in the introduction expends hundreds of millions of floating 
point operations to process each electronic orbital, each of which are represented in terms of only 
six thousand coefficients. 

Workers interested in electronic structure therefore have sought different methods. One ap- 
proach has been to take techniques from the wavelet literature such as the "non-standard" mul- 
tiply approach of Beylkin, Coifman and Rokhlin[BCR91], which has been applied to the solution 
of Poisson's equation in multiresolution bases with the processing of some additional coefficients 
but still leading to an efficient scheme[GI9£]. Workers also have developed new methods specifi- 
cally for physical calculations [LAE9J]. These new methods allow operators to be applied to the 
electronic wave functions directly in their "packed" representation without processing any addi- 
tional information and have been shown to be several times more efficient than the non-standard 
matrix approach in situations typical of electronic structure [LAE98|. The associated transforms 
have been shown to have the novel property that the process of (a) unpacking the physical fields 
at a number of points in space equal to only the number of packed coefficients, (b) coupling the 
physical fields in any local, non-linear fashion at these points, and (c) repacking the result always 
yields coefficients identical to what would be obtained with a fully unpacked function on a grid of 
arbitrarily fine resolution [LAE98]. With these latest advances, the field now stands poised to see 
the first applications of multiresolution analysis to large-scale electronic structure calculations. 



Ill Multiresolution Analysis of Electronic Structure 
III-A Multiresolution analysis and restriction 



Section IV reviews the mathematical structure of multiresolution analysis in detail. Here, we give 
a conceptual overview sufficient to discuss the use of multiresolution analysis in the calculation of 
electronic structure. 

Figure |] illustrates the concept of multiresolution analysis and its application to electronic 
structure. Throughout this work we shall use the variables Q, R and P to denote different levels 
of resolution and will take the coarsest and finest levels of resolution in a given calculation to 
be M and N, respectively. As [Mal8£, Mey90, Dau92, Chu92a. [5N96 | describe, a multiresolution 



analysis begins with a basis of coarse resolution Q = M, which consists of basis functions laid out 
on a regular grid across the region of interest (larger circles in the figure). The span of this set of 
functions, the vector space of all functions formed by their linear combinations, is denoted Vm- For 
reasons which will become clear below, the basis functions of this coarse space are referred to as 
the scaling functions. 

The next conceptual step is to increase the resolution of the basis by adding finer resolution 
functions at the points of a finer grid (smaller circles in the figure). The basis consisting of both the 
original coarse scaling functions and the new finer functions now spans Vm+i, a space of increased 
resolution Q = M + 1. The added functions are referred to as the detail functions or the wavelets. 

One may continue adding finer levels of detail functions to reach the final desired level of 
resolution, Q = N. We shall designate as Wq + \ the space spanned by the detail functions which 
bring the resolution from level Q to the next level Q + 1, so that Vq + ± = Vq © Wq + \. (Note 
that some authors prefer to designate the above detail space as Wq, rather than Wq + \.) Here, as 
throughout this work, by the addition of vector subspaces "©", we mean the space of all vectors 
which may be written as a sum of a pairs of vectors, one from each subspace. The space of all 
functions which can be described by functions on all of the scales included in the basis is thus 

Vn = V M © W M +i ® • • • ® W N . 
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Figure 4: Multiresolution analysis and the application of restriction to electronic structure: coarse 
grid (larger circles), detail points of finer grid (smaller circles), atomic nucleus (diamond), spheres 
of resolution (large circles centered on the nucleus), basis functions restricted from basis (empty 
circles), surviving basis functions (filled circles). 



Figure || shows the behavior of the expansion coefficients in such a multiresolution basis for the 
la state of the nitrogen molecule, which appeared in Figure |2| in the introduction. Figure || displays, 
on a logarithmic scale, the magnitude of the expansion coefficients for this state as a function of 
the two parameters characteristic of each basis function: location in space r and resolution Q. The 
three dimensional location r is projected onto the one-dimensional horizontal axis r as the distance 
from the center of the basis function to the nearest atomic nucleus. The scales Q of the basis 
functions are coded by different symbols. As evident in the figure, the separation of information 
into different length-scales afforded by the multiresolution analysis results in separate characteristic 
exponential decay envelopes with distance from the nuclei for the coefficients of each scale. These 
envelopes illustrate the fact that the finest scale coefficients need only be kept for basis functions 
in the immediate vicinity of the nuclei. This is precisely the behavior which makes multiresolution 
analysis so attractive for the calculation of electronic structure. 

The strategy introduced by |CAJL93[| to exploit this behavior is illustrated in Figure |||: about 
each nucleus we draw a set of successively inscribed spheres of appropriate radii for the scales 
M, M + 1, . . . , N, and we keep in the basis only those functions of a given scale which fall within 
the corresponding sphere. The grid points whose associated functions appear in the final basis 
set according to this prescription appear as filled circles in the figure. (Calculations with periodic 
boundary conditions generally include all functions on the coarsest scale.) We refer to this process 
of selecting grid points and their associated functions as restriction. For other applications, it is 
not always known a priori how to restrict the basis. For adaptive restriction approaches, the reader 
may wish to consult gPTgg , fJTWj , pK97| , |BK97| J. Tymczak and Wang |TW97| ] have also developed 
an adaptive restriction technique specifically for electronic structure calculations. 

Selecting the cutoff spheres so as to discard only coefficients below a given tolerance gives a 
systematic procedure for working with a dramatically reduced number of coefficients while main- 
taining a description equivalent, for any given tolerance, to the full basis Vjsr- In contrast to finite 
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Figure 5: Coefficients of the multiresolution analysis of the la orbital of N2: radii of restriction are 
denoted by vertical lines. 



element approaches, in multiresolution analysis there are no basis set derivative corrections in the 
Hellman-Feynman theorem because the basis functions of a multiresolution analysis remain fixed 
in space as the atoms move. The only effect of the motion of the atoms is to turn on or turn off 
basis functions whose coefficients are below the selected tolerance. As a result, the discontinuous 
effects from such on-off switching events are controllable and generally quite small. For example, 
in the calculations of [CAJL93], although the basis was chosen only to produce correct total ener- 
gies without regard to the calculation of forces, the small jumps in the forces calculated from the 
Hellman-Feynman theorem were only on the order of 1 meV/A. 



III-B Algebraic structure 

The main step in implementing a density functional calculation is to express each term of the LDA 
Lagrangian (|l|) in terms of the coefficients d a and c Qj j for the Hartree field 4>(r) and the Kohn-Sham 
orbitals {ipi(r)}, which appear in the expansions 

Mr) = $> a A(r) (4) 

a 

4>( r ) = ^dab^r), 

a 

where {b a (r)} is the basis set used in the calculation. Once the Lagrangian is represented in terms 
of these coefficients, the gradients of the Lagrangian with respect to d a and c a ^ may be calculated 
so as to locate the saddle point and thereby determine the orbitals {V'i(r)}, potential <f>(r) and total 
electronic energy C of the system. Several different approaches for expressing the Lagrangian have 
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been used in the application of multiresolution analysis to electronic structure. Here we present an 
overview of these approaches and the operators and transformations which they involve. Explicit 
details and formulae are given in Sees. III-D-QI-E. 

The simplest terms in the LDA Lagrangian are the electronic kinetic energy T and the Hartree 
field self-energy Vh-h, the first and final terms of (|lj). These are bilinear in the coefficients and may 
be evaluated exactly in terms of two-center integrals of the Laplacian operator between the basis 
functions. Because the matrix elements of the Laplacian operator are known explicitly for functions 
making up multiresolution analyses (Sec. V-C.2D , all implementations of multiresolution analysis to 
date use this simple two-center form to evaluate T and V g _ H [pAJL93| , |ACLT95| , |WC96| , |TW97[| . 

The electron-ion potential energy V e -i, the second term in Clda, is also bilinear in the expansion 
coefficients c a ^, but the matrix elements of the ionic potential needed to compute V e -i consist of 
three-center integrals, each involving two basis functions and one ion, and thus require special care. 
To date, V e -i has been handled through the introduction of a grid G of points p in real space. 
Two prescriptions for employing this grid exist. The first prescription (used in | ACLT95| ) follows 
the spirit of the energy functional and computes the electron-ion interaction as a functional of the 
electron density. This recipe first uses a forward transform to determine the values of the wave 
functions on the grid. From these, the single particle density n(r) = J2i f\' l l ; i( r )\ 2 is computed on the 
grid and the corresponding expansion coefficients are determined through an inverse transformation. 
From the resulting expansion coefficients, the final step is to compute the total potential energy in 
terms of known overlaps between the basis functions and the ionic potential. This approach has the 
advantage that the grid need only be able to resolve the electron density n(r) and not necessarily 
the potential Vj on (r), which may vary much more rapidly. We refer to this approach below as the 
energy functional prescription. 

The second prescription (used in [WC96, TW97| ]) follows the spirit of the effective Schrodinger 
equation for the Kohn-Sham orbitals and applies the electron-ion interaction as a diagonal operator 
in real space. This recipe also begins with a forward transform to evaluate the wave functions ipi{r) 
on the points of the grid. It then applies the potential operator at each point p of the grid G to 
produce V- lon {p)ij)i{p) and inverse transforms the result to coefficient space. Finally, the overlap 
of ip*(r) and V^ Qn (r)ipi(r) is computed using the known overlaps between the basis functions. We 
refer to this approach below as the operator prescription. 

The fourth term of @) describes the coupling between the electrons and the Hartree field, V e -u- 
This term is cubic in the expansion coefficients and involves three-center integrals. In principle, 
a treatment in terms of direct three-point interactions is possible using analytic results for the 
integrals of triple products of scaling and detail functions developed in BK97]. To date, this direct 
route has not been pursued in the calculation of electronic structure. Instead, one may exploit the 
fact that the coupling V e -H has the same structure as the electron-ion coupling and compute it in 
the same manner, either as a functional of the electron density or as the result of the application 
of 4>(r) as a local operator. 

The final remaining term, the third term in (|l|), gives the local density approximation to the 
total exchange-correlation energy of the system, E xc . This is known only in terms of the non- 
algebraic function e xc (n), which is tabulated in [PZ81], and thus cannot be evaluated in terms 
integrals of products of basis functions using the formalism of [BK97]. For this term, there is 
no choice but to evaluate the electron density on a grid G and evaluate the tabulated function 
£xc( n ) on a point by point basis in real space. Once this is done, there are a variety of choices for 
how to use the exchange-correlation field to determine the total exchange-correlation energy E xc . 
One could proceed directly and evaluate the exchange-correlation energy density per unit volume 
e xc( n (p)) ■ n (p) a t each point and then integrate the result numerically, which would introduce an 
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additional class of approximation in the evaluation of the Lagrangian. What has been done, instead, 
is to follow the same prescriptions as used in computing the electron-ion interaction, following either 
the energy functional prescription and computing E xc as the integral of the product of n(r) and 
£xc(^-(^*))[ |ACLT95 ] or following the operator prescription and computing E xc as the expectation of 
e xc (n(r)) as an operator acting on the Kohn-Sham orbitals[ WC9C , TW97]. 



III-C Matrix language 

To discuss the various approaches available for evaluating the LDA Lagrangian in terms of the 
expansion coefficients for the physical fields, we now introduce a matrix language for electronic 
structure methods which are based on single-particle orbitals. In analogy with Dirac's bra-ket lan- 
guage, this matrix language is completely explicit but keeps the expressions for physical quantities 
independent of the details of the underlying basis set. This allows us to discuss the various strate- 
gies for applying multiresolution analysis to electronic structure without obscuring the discussion 
with irrelevant details. This new language is useful in its own right for both formal manipulations 
and the development of highly portable, efficient software. Plane wave calculations performed with 
software developed by translating this language directly into the C++ include | IBA98| , |CIBA98|1 . 



For simplicity, we consider below cases where the Fermi occupations / are constant and the 
sampling of the Brillouin zone is carried out at the T point. More general cases including variable 



fillings and multiple k-points may be worked out with some additional complexity [ APJ92 1 . For 
systems sufficiently large that the Brillouin zone may be sampled at the T point alone, the formulae 
below may be used directly. 



III-C. 1 Fundamental basis-dependent operations 

We begin by introducing two operators, diag and Diag . The operator diag converts a matrix to 
a column vector containing the elements along the diagonal of the matrix, and the operator Diag 
converts a column vector to a diagonal matrix with the components of the vector placed along the 
diagonal. In terms of components, for a matrix M and a vector v, these operators are 

(diagM) Q = M aa (5) 
(Di&gv) a/ 3 = v a 5 a/ 3, 

respectively, where 5 a p is the Kronecker 5. Note that while diag Diagu = v for any vector v, 
Diag diag M = M if and only if the matrix M is diagonal. Two identities involving these operators 
which we shall use freely are 

(diagM)^ = Tr (M+Diagv) , (6) 
w f diagM = Tr ((Diagv^Af) . 

Next, it is useful to regard the Hartree field expansion coefficients, the d a from (Q), as the 
components of a column vector d and the electronic wave function coefficients, the c a ^ from (||), as 
the elements of a matrix C, each of whose columns contains the expansion for a single electronic 
orbital, C a i = c a ^ For formal manipulations, it is convenient to define also 

P = fCC\ (7) 

the representation of the single particle density matrix in the space of basis functions {b a (r)}. 
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O a(3 = J d\b* a {r)bp{r) 
L aP = f cPrb* a {r)V 2 bp{r). 



As Sec. [II-E describes, evaluation of the Lagrangian (|l]) in terms of the expansion coefficients 
contained in d and C requires knowledge both of the overlaps of the basis functions among them- 
selves and of their matrix elements through the Laplacian operator, 

(8) 
(9) 

The first two non-trivial basis-dependent operations which require careful implementation are mul- 
tiplication by the inner product matrix O and by the Laplacian matrix L, to which we refer 
below as application of the overlap operator and application of the Laplacian operator, respectively. 
Section V-C.2| discusses how to compute the required matrix elements, and Sec. VI describes 
efficient techniques for applying these operators. Note that in the special case of orthonormal 
bases [ WC9G| , TW97], we have simply O = I, where I is the identity matrix. 

The forward transform operation described in the previous section converts the expansion coef- 
ficients of a function into the values of the function on the points p of the grid G in real space. This 
operation simply amounts to multiplication of a column vector containing the expansion coefficients 
by the matrix 

Z pa = b a (p), (10) 
whose aft 1 column consists of the values of the afi 1 basis functions at all of the points p of the grid. 



In the case where different bases are used for the wave functions and for the Hartree field | WC96 , 



TW97], the columns of X consist of two subsets, one for each of the two basis sets. 



As described in Sec. III-B, it is at times necessary to find the expansion coefficients for a 



function from its values on the grid G. We denote this linear inverse transform operator as J . In 
implementations where the number of grid points equals the number of basis functions [ |ACLT95 , 
LAE98j , Ycs97|, the natural choice is to take J = Z^ 1 . However, we maintain the distinction 
between J and Z~ l because in implementations where a full uniform grid G of sampling points is 



used |WC96| , |TW97|| , there are generally more grid points p than basis functions b Q (r), and Z and 
J cannot be inverses. In the case where more that one basis set is used[WC9£, TW97H , formally, J 
computes the inverse transform separately for each basis set. In practice, however, as seen below, 
only one inverse generally need be computed because J will be needed only on one basis set. 

Finally, as we shall see in Sees. 1II-D| and QI-E , two further transforms appear in the calculation 
of the gradients of the Lagrangian. These two conjugate transforms represent multiplication by the 
Hermitian conjugates Z^ and J* of the standard transforms. The reader should bear in mind that 
the relation Z^ = Z^ 1 = J for the discrete Fourier transform is quite special. Because orthogonality 
with respect to integration does not ensure orthogonality with respect to a discrete sampling of 
the functions the above relation for the discrete Fourier transform is not generally true, even for 
orthonormal bases. In particular, Z^ ^ Z~ l for the multiresolution bases of Daubechies wavelets 



used in [|WC9q , p?W97fl . 

In summary, six non-trivial basis-set dependent operations are needed in the calculation of 
electronic structure: the application of two operators (the Laplacian and overlap operator) and 
four transforms (forward, inverse, and the conjugates to each), which we denote as L, O, I, J ', Z^ 
and j\ respectively. 



III-C.2 Identities 

Although the action of the six operations L, O, Z, J, Z^ and depend on the specific choice of 
basis, they obey several important identities of which we shall make use. In addition to their use in 
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formal manipulations, these identities provide a useful practical tool to verify the implementation 
of the various operators. 

The most important among these identities involve the constant function. To represent the 
constant function on the grid, we introduce 1, a column vector containing unity as each entry on 
the grid: 

h = l. (11) 

Plane waves, finite element bases, and proper multiresolution analyses all have the property of 
being able to represent this function exactly. For these bases, we have that for all points r, 
^l a {J^)aba{r) = 1. Evaluating this on the points p of the grid G yields that, in particular, 

1J1 = 1. (12) 

For other sufficiently descriptive bases, such as Gaussian bases, this and the relations below should 
hold at least approximately in the regions described by the basis. 

There is a close relationship between the inner product matrix O and the integrals of the basis 
functions. The vector 

s = OJ1, (13) 
is a column vector containing the integrals of each of the basis functions: 

s a = (OJl) a 

d 3 r b* a (r) \Y,(Jl)pbp(r) 



= Jd\b* a (r). 



Thus, if g is a vector of expansion coefficients, the integral of the function represented by g is 

f d 3 rg(r) = f d 3 r ^g a b a {r) = s^g. (14) 

From this, we may also derive the normalization condition 

s ] Jl = [ld 3 r = Q, (15) 



where is the volume in which the calculation is carried out. 

In solving Poisson's equation for the Hartree potential, care always must be taken with the null 
space of the operator L. Integrating the identity V 2 1 = against the complex conjugate of each 
basis function gives 

LJ1 = 0. (16) 

In periodic systems, where the only solution to Laplace's equation is the constant function, the 
entire null space of the operator L consists of the vector J\. Below, we use this to determine the 
precise value of the compensating average density no needed to avoid divergences in the Lagrangian 
in periodic calculations. 

Finally, although there is no a priori relationship between and J", an approximate relation 
exists when the grid G is uniform and of high resolution. Under these conditions, we have 

OJ w ul\ (17) 
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where uj is the volume per grid point. To see this, consider two arbitrary functions g(r) and h(r), 
where g(r) is represented by the vector g of its expansion coefficients and h{r) is represented by the 
vector h of its values on the points p of the grid G. There are then two ways of approximating the 
integral / d 3 r g(r)*h(r), both of which should give nearly the same result. First, one could inverse 
transform h to find appropriate expansion coefficients and then evaluate the overlap in coefficient 
space, giving the result g^OJh. Alternately, one could determine the values of g on the grid by 
the forward transform and then approximate the integral as the sum u(Igyh. Equating these two 
expressions for general vectors g and h leads to ([l7|). 



III-D Lagrangian and energy functionals 

As described above, the cornerstone of all density functional calculations which use basis set ex- 
pansions is the explicit expression of the LDA Lagrangian in terms of the coefficients d and 
C of the expansions (0). The language and operators defined in the previous section were created 
specifically so that the expression of the Lagrangian is identical for most basis sets, including plane 



wave, Gaussian orbital, and multiresolution bases. As described in Sec. III-B two common strate- 
gies exist for expressing the Lagrangian in terms of the expansion coefficients in multiresolution 
analyses, the energy functional prescription and the operator prescription. We now give explicit 
expressions for these two prescriptions. 



III-D. 1 Energy functional prescription 

III-D. 1(a) Lagrangian The only terms in the Lagrangian (|l|) which cannot be written directly 
as a functional of the electron density are the electronic kinetic energy and the Hartree field self- 
energy. These two terms are simple bilinear forms in the expansion coefficients and are best 
evaluated exactly by the direct substitution of the expansions @ into (pi), giving 

T = -lT,fT,^ L ^ = -\^f^LC=Tr((-h)p) (18) 

i a/3 

V H -h = -^-y2d* a L al3 dp = ^Ld, (19) 

S7T T 87T 

ap 

where we have detailed the explicit conversion to matrix language and made use of cyclic property 
of the trace to give a representation of T in terms of the density matrix P defined in (fjj). 

The remaining terms in the Lagrangian all may be written directly in terms of the electron 
density. In our matrix representation, each column of the matrix IC contains the values of one of 
the electronic orbitals when evaluated on the grid. The electron density at each point p of the grid 
is thus 

n(p) = J2f( IC U IC )pi 

i 

= (/(TC)(JC)t) . 

\ / pp 

From this, we may gather the real space charge density into the column vector n as 

n = diag (jPJ t ) . (20) 
Finally, the expansion coefficients for the function n(r) are just Jn. 
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Using this last result to evaluate V e -i, we find 



d 3 r Vj on (r)n*(r) 
(Jn^v 

Tr (rf(Dmgjiv)lP) , 



(21) 



where the vector 



<Prb* a (r)V ion (r) 



(22) 



represents the ionic potential as its overlap with each function of the basis. (In conjunction with 
the appearance of no and the use of the Ewald summation, the zero wave-vector component of 
V[ on should be subtracted when evaluating these overlaps in calculations with periodic boundary 
conditions.) Although n(r) is real, we introduced n*(r) above as a formal device to reduce the 
number of complex conjugations appearing in the final expression. The conversion to density- 
matrix form in the last line of ( ^1~1) may be carried out using the identities @. 

The evaluation of the next term in the Lagrangian, the total exchange correlation energy E xc , 
is similar. Given access to the values of n on the grid points, it is a simple matter to also evaluate 
the exchange-correlation energy per particle e xc (n(p)) at those points. Collecting these into the 
column vector e xc (n), inverse transforming both this and the charge density and taking the overlap 
in coefficient space gives the final result, 



E xc = J d 3 re xc (n{r))n*(r) 
= (Jn)^0{Je xc (n)), 



(23) 



where again we introduce the conjugation of n to reduce the number of complex conjugations 
appearing in the final result. E xc may also be converted to density-matrix form using @. 
The final term remaining in (||), V e -H, also has the overlap form 



e-H 



d 3 r(j)(r) (n(r) — no)* 



(24) 



-8fc \{J(n - n l)) t Od] , 



where 1 is as defined in ([n]). Again, the conjugation of the density term has no effect but to yield 
the simplest final expression. 

To determine the proper choice of no for periodic supercell calculations, we note that the Hartree 
self-energy term Vh-h has no contribution from the projection of d in the null space of L, which 
Sec. III-D| showed lies along the direction J\. Thus, in order for a saddle point to exist for the 
Lagrangian (|), there can be no coupling of this component of d to the electron density in (|24|). 
Hence, we must have {Jin — nol))^0 • = 0. With the identities of Sec. [III-C.2; , this means 
that no = \Jn) /Q, in accord with our interpretation of no as the integral of n(r) divided by the 
volume of the supercell. With this result, V e -H is 



Ve, 



H 



-5? 



J ] {0- S -^-)d 



n 



n 



(25) 
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Putting the preceding results together, the final expression for the LDA Lagrangian under the 
energy functional prescription is 



£lda(C, d) 



Tr [f C \--L)C)+(Jn^v 



+{Jri) ] OJe xc {n) 



ss 1 



n 



8-7T 



Tr (j-^L)p) + Tr Diag J^v)IP 
+ Tr (jt ( Diag J j OJe xc (n) ) IP 



TrX 1 ^ 



Dia gl 7t (0 -_)d 



lP+—d^Ld, 
8?r 



(26) 



(27) 



where n is defined as in (|2(j). Here we have given two forms. The form (^) is more efficient for 
evaluating the value of the Lagrangian in terms of the coefficients C and d. The form ( |27| ) , which 
is expressed almost entirely in terms of the density matrix, is most useful in formal manipulations 
which determine the gradients of the Lagrangian with respect to the electronic coefficients. 



III-D.l(b) Poisson's equation To determine the expansion coefficients of the Hartree field (p, 
one sets to zero the variation of the Lagrangian with respect to all possible infinitesimal variations 
5d in the Hartree field coefficients. Taking the real part of ([^) explicitly, this variation is 



5£lda 



1 



•(")M5(°-^» + s m ; 



i 



(28) 
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t 



SS ' 1 

- 7T )Jn + —Ld\ Sd. 

Si 07T / 



The gradient with respect to the real and imaginary parts of d are thus the real and imaginary 
components of 



ddC-LDA 



SS' 



(29) 



respectively. Setting this to zero, we arrive at the basis independent representation of Poisson's 
equation, 



ss 



Ld = 4ir(0 - —)Jn. 



(30) 



The positive sign of the right side of this equation reflects the negative unit charge of the electron. 

Care must be taken in solving this equation because, as described above, L has a null space 
along the direction J A.. However, by the construction of (p5|), the right-hand side of the equation 
has no projection in this space. Thus, in the solution 



col 

d = 4irL- 1 (0- —)Jn, 



(31) 



L is understood to mean inversion of the linear operator L in the sub-space orthogonal to the 
null space and leaving zero projection along the null-space in the result. 
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III-D.l(c) Energy functional Substituting the explicit result for d into the Lagrangian ([]]) 
gives the following explicit expression for the LDA energy functional Elda(C), 

E LDA (C) = Tr(ftf(~L)C\ + (Jn?v (32) 
+{Jri) ] OJe xc {n) 

+\n^j\0 - ^)(-4^L~ 1 )(0 - a A) Jn . 
Note that Hermitian symmetry has ensured that the Hartree term is explicitly real. 



III-D.2 Operator prescription 

III-D.2(a) Lagrangian The alternative prescription for expressing the Lagrangian, followed 
in [ WC96| , [TW97 1 , is to view it as the sum of the expectations of an energy operator C among the 



orbitals, Clda = Z)i / / d 3 ril>l£ipi. The only term which cannot be taken to involve involve the 
electrons in this way is the Hartree self-energy term, which is best represented in the exact form 
(|l^) . Note that we have already expressed the electronic kinetic energy in its operator form in (|i~8|). 

The remaining contributions to e are local operators in real space and thus all take on the 
same form as does the electron-ion interaction V e —i- The technique which Wei and Chou| WC96f 



introduced to treat V e _j is to first compute ZC, the values of the wave functions on the grid points 
and then multiply the result by a diagonal matrix containing the values of the potential at the 
grid points. This results in (Diagu)ZC, where v is the vector containing the values of the ionic 
potential at the grid points, 

(v) p = V ion (p), 

where the tilde on v distinguishes this vector of values from the vector v of overlaps defined in ( P2|) . 
To compute the energy, this result is transformed into expansion coefficients by the operation of J 
and the overlaps with the ip*(r) are computed using the known overlaps between basis functions. 
The total of the resulting potential energy among all of the orbitals is 

V e . % = E// d3r ^(0(^ion(0^(r)) (33) 

i 

= Tr ftf • O ■ J ( ( Diag v)lC) 
= Tr ftfVC = Tic VP, 

where V = OJ{ Diagu)T. A potential difficulty with this approach is that the inherently asymmet- 
ric roles played by if}* and i/j have resulted in a potential energy operator V which is not Hermitian 
and thus may lead to complex energies. Taking the Hermitian part of V corresponds to taking just 
the real part of (^) and gives the proper form for the electron-ion energy. After some manipulation 
and using the fact that v is real, we have the following equivalent forms for V e -i 

V e -i = KTr(O l 7Diag{)T-P) (34) 
= «t»diag (IPO J) 
= uiv^h, 

where we have defined 

n = $tdiag(lPOJ)/u, 
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as an effective charge density, and u as the volume per grid point. With these definitions, the 
electron-ion interaction takes on precisely the appearance of a numerical integration in real space 
of the product of potential and the effective electron density. Indeed, from ( |i~7| ) and (|20|), as long 
as a fine, uniform grid is used, we have 3ft diag (IPO J) ~ urn, so that h approximates the physical 
electron density. As we shall see, the effective density h replaces the physical density in all energy 
expressions of the operator prescription. 

The remaining terms from (jl|) have the same structure as V e -i and thus may be evaluated in 
the same way to yield 



E, 



v; 



e-H 



-un[{Id)^{n - n l)]. 

To determine no, we again ensure that there be no coupling between the electron density and the 
projection of d along the null space of L: (I ■ J"L)\n — uqI) = 0. This simplifies to 

1+n 
no = TfT , 

the mean value of h among the points of the grid. 



(35) 



We shall find in Sec. III-E.2 that the the exchange-correlation energy should also be evaluated 
using the effective density n, so that the final expressions for the Lagrangian in the operator 
prescription are 



£lda(C, d) 



Tr (fC\-~L)C \ +ujtfn 
+w(e xc (n)) t n 



(36) 




+ ^-Shd 



Tr l^(-^L)Pj +5RTr (OJ{ Diag£)XP) 
+5RTr (O J Diag (e xc (h))lP) 

-KTr ( OJ I Diag 5? 



(37) 



ii-^id^xpyi-SLd. 



III-D.2(b) Poisson's equation Following the same procedure as d28|) to compute the variation 



of (36) with differential changes in the Hartree field coefficients 5d gives, in the operator prescription, 

d d C LDA = [I - ^7 ) h + ^-Ld. (38) 



Thus, 



Ld 




(39) 



is the representation of Poisson's equation in the orbital approach. Note that the effective electron 
density n, rather than the actual density n is the source term. 

In solving this equation for d, the inversion of L _1 is again understood to take place in the 
space orthogonal to the null space of L, giving the result, 

lit 
ltl 



d = AttuL- 1 !^ [i - -r- 



n. 



(40) 
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III-D.2(c) Energy functional With the result (Sc|), the expression for the LDA energy func- 
tional Elda(C) in the operator prescription is 



Elda{C) 



2 

oj(e xc (h))^n 



(41) 



urn. 



III-E Kohn-Sham equations 



In the preceding sections, we computed the gradient of the LDA Lagrangian fl2q,36) with respect 
to the Hartree field coefficients to derive the Poisson equation. In this section, we consider the 
gradient with respect to the electronic coefficients to determine the effective Schrodinger equation 
for the electronic orbitals. 

Before proceeding, we first note that as a consequence of the Hellman-Feynman theorem, the 
expression for the gradient with respect to the electronic coefficients of the LDA Lagrangian may 
also be used for the gradient of the LDA energy functional, so long as one substitutes the appropriate 
expression ([3l] or |40| ) for the Hartree potential coefficients in terms of the electronic coefficients C. 
This follows from the fact that the solutions ({31,40) ensure ddJ~-LDA(C,d(C)) = and thus 



^E LDA {C) 



—C LDA (C,d(C)) 

d c c LD A(c, d(c)) + d d c LDA (c, d{c)) ■ d c d{c) 
d c c LDA {c,d{c)). 



(42) 



It thus suffices for us to determine the gradient of the Lagrangian with respect to C . 

This is done most efficiently by expressing all contributions to the variation of C-lda in the form 
Tr M5P where M is a Hermitian matrix and 5P is the variation in the density matrix defined in 
(|7|). In terms of the variations in C, such a variation takes the form 



Tr (M 5P) = Tr ((/A/C) f 5&j + Tr Utf (fMC) 



(43) 



The contribution from such a variation to the total gradient with respect to the real and imaginary 
parts of of C are therefore the real and imaginary components of 2/MC, respectively. 

The final consideration when varying the wave function coefficients is the orthonormality con- 
straint ©, which in our matrix language becomes 



C^OC = I. 



(44) 



The analytically continued functional approach [ APJ92 | deals with these constraints by introducing 
a set of orbital expansion coefficients Y which are unconstrained and may have any set of overlaps, 



U = Y^OY. 



(45) 



The coefficients C are then defined as dependent variables through the mapping C = YU^ 1 ^ 2 , 
which ensures that the constraints (El) are always satisfied automatically, as easily verified by 
direct substitution. 
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In terms of the independent variables Y, the density matrix appearing in the functionals ( f27| , [3^ ) 

is 

p = fCC^ = fYU- x Y^. (46) 
After some manipulation, the variations with respect to the Y take the form 

TrM5P = Tr M8(fYU" x Y^) (47) 
= Tr (/PMyC/^ 1 ) f 5Y + Tr 5Y^ (/PMyf/^ 1 ) , 

where we have used the relation 5{U~ l ) = —U~~ 1 (5U)U~ 1 and defined 

P= (I-OYU-W) = (I-OCCft), (48) 

an idempotent projection operator onto the subspace spanned by the unoccupied wave functions. 
From these considerations, we have that the contributions from such a variation to the gradient of 
£>lda with respect to the real and imaginary parts of the unconstrained variables Y are the real 
and imaginary parts of 2fPMYU~ 1 , respectively. 



III-E.l Energy functional prescription 

All but one of the dependencies of ( p7|) with the wave functions are already explicitly in the form 
Tr MP, so that their contribution to the gradient may be evaluated immediately according to 
). The one remaining dependency of (p7|) on the electrons is through the density dependence 



of 

£xc{ n )- This variation may also be cast into the form of ( |43| , |47 |). To do so, we consider the effect 
of this variation on E xc when expressed in the form (^) . Using the definition of n (^) and one of 
the identities @, this leads to 



rJj ] OJ5e xc {n) 



n 



] J ] OJ ({V)\&ge' xc {n))8n) 



Tr 



The total variation of (| 



jtDiag (( Diag e' xc (n))J^OJn^ 15P 
as the wave functions vary is thus, 
5C LDA = Tr (H sp 5P) , 



where 



H sp = -ix + TtDiag {jW J^OJe xc (n) (49) 
+ (Diag e' xc {n))J^OJn 

is the energy functional prescription for the representation of the single-particle Kohn-Sham Hamil- 
tonian. Using (43,47), the gradients which we seek are therefore 



dcC LDA (C, d) 
d Y C LDA (Y,d) 



2fH sp C 
2fPH sp YU-\ 



(50) 



where U and P are defined in ( fi"5| ) and (|48|), respectively. 
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III-E.2 Operator prescription 



The results for the orbital approach are directly analogous. All of the dependencies of the La- 
grangian in the operator prescription on the wave functions, except for that of e xc , are already 
explicitly in the form TV MP. If we take e xc as a function of n, then this variation of E xc (second 
line of ©) is 



w(<5e a;c (n)) t n 



o;(Diage^. c • 5h)^h 

uj ((Diag e^Jn) 1 " 5n 

U Tr (OJ{ Diag ( Diag e' xc ■ h))T8P) 



This combines with the explicit dependency of the exchange-correlation term in (37) on P to 
produce the total exchange-correlation energy dependency in the form 



6E a 



8(ue xc (n)^ h) 

KTr (OJ(Biagv xc )15P), 



where v xc is the usual exchange-correlation potential operator, 

v xc = dn(e xc (h) n) = e' xc {n) n + e xc (h). 



(51) 



Note that, as referred to in Sec. II 1 -D.2| this form is obtained only when the exchange correlation 
energy per particle e xc is evaluated as a function of the effective density h, not as a function of the 
physical density n. 

Combining this result for E xc with the remaining dependencies in the operator prescription for 
the Lagrangian ( |37| ) gives the total variation in the form SClda = TrH sp 5P, where 



H 



sp 



1 r 1 

■-L+- 



OJ Diag lv+ [(Diage^n)) n + e xc {fi)] - U 





lit 








| J + h.c. 



, (52) 



and "h.c." stands for the Hermitian conjugate of the entire matrix product to the immediate left. 
Again, once given this expression for H sp , the expressions (50) gives the gradients with respect to 
the wave functions. 



III-F Solution techniques 

Locating the saddle point of the LDA Lagrangian (|l]) corresponds to the simultaneous, self- 
consistent solution of the Poisson and Schrodinger equations. Although the Lagrangian formulation 
allows for a direct search for the saddle point, most calculations to date have solved Poisson's equa- 
tion explicitly at each iteration and then followed the gradients in ( |50| ) to solve the Schrodinger 
equations using standard electronic structure approaches. The approaches applied to solve the 
Schrodinger equations include conjugate gradient minimization with respect to y[ACLT95] and 
iterative matrix diagonalization[ WC96] and Car-Parrinello simulated annealing [TW97] to find C. 
The solution of Poisson's equation in multiresolution analyses, on the other hand, has garnered 
special attention. 

To solve the Poisson 



Wei and Chou[WC96] introduced the innovation of using a separate 
Fourier representation for the operator L in the self-energy term of the Hartree field. The inter- 
pretation of (|40j ) when using such a representation is the following procedure. First, subtract the 
average value from the effective density n, then perform a discrete Fourier transform of the result 
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Convergence of Preconditioned Conjugate Gradients for N2 




40 60 80 

Number of iterations -> 



120 



Figure 6: Convergence of iterative conjugate- gradient solver used in [ ACLT95[ | (results from 
| LAE98| J). 



into Fourier space, divide by the familiar factor of negative squares of wave vectors — q 2 (ignoring 
the q = term), and scale by Attlo to obtain the final result. The disadvantage of this procedure is 
that one must store and process data associated with an unrestricted grid of points at the highest 
resolution. 

When working only with the points of a restricted grid, no such direct solution is known. 



Instead, iterative solutions have been used. In [ACLT95], a conjugate gradient solver with simple 
diagonal preconditioning was used to solve Poisson's equation in a semicardinal multiresolution 
analysis of the type described in Section V-B.3 . Figures || and show the results of a study of 
the performance of this algorithm [LAE9£]. Figure shows the root mean square magnitude of the 



residual as a function of iteration as this algorithm is applied to the nitrogen molecule in a basis 
with seven levels of refinement. After an initial phase of about twenty iterations, the convergence 
of the solution of Poisson's equation becomes nearly perfectly exponential and reduces the residual 
by over ten orders of magnitude in one hundred iterations, very good performance for a system 
consisting of 14,000 degrees of freedom in which the Laplacian operator has a nominal condition 
number of 2 x 10 6 . 

The slope of the exponential portion of the convergence curve in Figure || corresponds to a 
reduction in the error at each iteration by 25%. Defining the effective condition number of the 
procedure as the inverse of this fractional improvement gives an effective condition number c ~ 4. 
Figure explores the behavior of this effective condition number as a function of the number of 
refinement levels k. Beyond about five refinement levels, the effective condition number remains 
essentially constant, implying that a constant number of iterations suffices to produce a result of 



a given accuracy. With the methods described in Sec. VI-B, the computational work involved in 
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Performance of Preconditioned Conjugate Gradients for N2 
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Figure 7: Effective condition number of iterative conjugate- gradient solver used in [ACLT95] as a 
function of resolution. (Results from | LAE98| .) 



each iteration is linear in the number of coefficients n r in the restricted multiresolution analysis. 
For these problems, this simple preconditioned conjugate gradient approach therefore produces the 
solution to Poisson's equation with 0{n r ) operations. 



Working with special linear combinations of the basis functions from [LAE98] known as lifted 



wavelets [Swe96], Goedecker and Ivanov have also used preconditioned conjugate gradients to solve 
Poisson's equation and report similar convergence rates and an 0{n r ) solution as well[GI98]. 

Finally, multigrid techniques have been combined with multiresolution analysis [ RJZ94 7 fycs97 ] . 
This approach has also been shown to produce 0{n) solutions of Poisson's equation in tests carried 



out in one dimension | Yes9 7 , YA9S 



IV Theory of Multiresolution Analysis 

In this section we review multiresolution analysis, the conceptual basis of wavelet theory. 
Our presentation below departs somewhat in notation and perspective from the traditional 
discussion [pau92] , Chu92a, [5N96 | to provide a presentation more suited to physical calculations 
in multiple dimensions. 

The basic idea behind multiresolution analysis is to provide a framework for systematically 
increasing the resolution of a basis in such a way as to always maintain a uniform description of 
space. Each stage of a multiresolution analysis doubles the resolution of the basis in the precise 
mathematical sense described in Sees. 1V-A and |IV-B . In order for this to be possible, the basis 
functions at the coarsest scale must be chosen from among a special class of functions known as 
scaling functions, whose explicit construction we describe in Sec. [V-Q . To analyze a function into 
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Figure 8: Bases involved in two-scale decomposition: coarse resolution space (Vo), fine resolution 
space (V\), detail space (Wi). 



contributions on separate scales and thus make possible the efficient representation of electronic 
structure described in Sec. Ill , a second type of basis function, known as either as a wavelet or a 



detail function, is needed to carry the finer scale information. For a given choice of scaling function 
for the coarsest scale, the use of only certain detail functions results in the systematic doubling 
of resolution at each stage in the multiresolution analysis. Sections IV- D and IV- E discuss this 
doubling process, known as two-scale decomposition, and specify the allowable classes of detail 
functions. The last paragraph of Section IV- E sketches an alternate, more basis-set oriented, 
viewpoint on two-scale decomposition which the reader may find illuminating. Once suitable scaling 
and detail functions are decided upon, including them into the multilevel framework described in 
Sec. IV- B| completes the construction of the multiresolution analysis. Our review of multiresolution 
analysis concludes in Section IV- F with the introduction a matrix representation which is useful 
for developing and describing techniques for electronic structure and other physical applications. 



IV-A Bases of successive resolution 

The top panel of Figure || illustrates a typical coarse resolution basis. In a multiresolution analysis, 
this starting basis consists of a set of d— dimensional functions <j>{x — n) translated to be centered 
on the points n of Z d , the cubic lattice of integer spacing in d— dimensions. The vector space of 
functions which may be expressed as linear combinations of these basis functions is denoted Vo- All 
functions F{x) in Vq thus may be expanded with expansion coefficients F n as 

F{x) = £ F n cf>{x-n). (53) 

neZ d 

Starting from Vq, one approach to produce a set of bases of successively higher resolution is to 
reduce the scale of space by increasing powers of two, compressing both the lattice on which the 
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functions are centered and the basis functions themselves. This process defines a sequence of lattices 
Cq,C\, . . . ,Cn, ■ ■ ■, associated with a sequence of bases of increasing resolution, Vq, Vi, . . . , V/v, . . ., 



c « s <54) 

Vq = span {</> (2 Q (x ~P))} peC = s P an {<P{^ Q x - n )} neZd ■ 

Note that the two representations for the basis functions of Vq given in the second line of (^) are 
equivalent. Below, we shall use whichever form is most suitable. The top and bottom panels of 
Figure § serve to illustrate the effect of going from Vq to V\. 

As we will see in Sees. |IV-C| — |IV-E| , it is possible under quite general circumstances to find 
additional basis functions, referred to as detail functions or wavelets, which when combined with 
the basis functions for the space Vq span the space of next higher resolution, Vq + \. The advantage 
of this approach for increasing the resolution is that it clearly separates the description into different 
length scales. The original coarse basis functions will continue to carry information about coarse 
scale behavior, and the new functions carry the higher frequency details. Accordingly, the space 
spanned by the detail functions is referred to as the detail space Wq+\. The center panel of Figure 
H illustrates the appearance of such a detail space. 

All functions in Vq+i may be decomposed into the sum of one function from the coarser space 
Vq and another from the detail space VFq+i, which means 

Vq(BWq +1 =Vq +1 , (55) 

where the addition is in this sense of combining vector spaces. We shall refer to this representation 
of the functions in Vq+i as the sum of a function in Vq and a function in PFq+i as the two-scale 
representation or the two-scale decomposition. We refer to the alternate representation of the same 
space functions in terms of the compressed scaling functions of scale Q + 1 as the single-scale 
representation. 

Clearly, for the single- and two-scale representations to be equivalent, the spaces Vq + \ and 
Vq © VFq+i must be of the same dimension. The basis spanning the finer space contains one 
function for each point in Cq+i, and the basis spanning the coarser space has one function for each 
point in Cq. Therefore, the basis spanning Wq + \ must consist of one function for each point in 

Dq + i = C Q+ i - Cq. 

In Figure |8| for example, Cq consists of all integer points, C\ contains all integer and odd half-integer 
points, and D\ consists of only the odd half-integer points. Figure ^ illustrates the appearance of 
these sets of points in d = 2 dimensions. 

These various periodic sets of points appearing in the discussion of multiresolution analysis may 
be described simply in the language of crystallography. The set -Dq+i contains the "decoration" 
points which when added to the lattice Cq produce the lattice Cq+i. Dq + \ therefore is a crystal 
of points with a crystalline basis containing 2 d — 1 points and with Cq as the underlying Bravais 
lattice. When dealing with a finite number of levels of resolution, Cm C Cm+i C . . . C CV, it 
proves convenient to have a notation for the set of points beyond a given level Q, 

B q = C n -Cq = D q+1 U...UD n . (56) 
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IV-B Multiresolution analysis 

Given suitable detail spaces Wq, iterative application of two-scale decomposition can augment a 
space Vm of fixed coarse resolution to a space Vn of arbitrarily high resolution, 

Vn = (v © Wi © . . . © W" M ) © w M +i © • • • © w N 

= V M © Wm+i © ■ • • © W N , 



where M < N in accord with the convention introduced in Sec. III-A. Correspondingly, all physical 



functions may be described by a sequence of evermore detailed components, 

f(x) = F M (x) + G M +i(x) + ... + G N (x) (57) 

where Fm(x) is a function in Vm and the Gn(x) are functions from the detail spaces Wq. We refer 
to the decomposition of a function in this form as the multiresolution analysis of the function. 



The analysis of functions in the form (57) is what gives rise to the great advantages of using 
multiresolution analysis. As illustrated previously in Figures Q and ||, the detail component func- 
tions Gq(x) are physically relevant only within the tiny volumes of the atomic cores and thus may 
be described by a dramatically restricted subset of the basis. 



IV-C Scaling functions 

For the condition ( |55| ) to apply, each basis function of the coarser basis Vq, including in particular 
the function c/)(2 < ^y) centered at origin, must be among the functions in the space Vq+i. Thus, 
there must exist some set of coefficients c n for which 

<K2 Q y) = E c n <X2 Q+1 y-n). 

neZ d 
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This condition is universal and independent of the scale Q as may be seen by the substitution 
x = 2®y, which re-scales the condition to scale Q = 0: 



<j>(x) =Y J Cn4>{2x-n). (58) 



Condition (|5§|) is known as the two-scale relation. Those functions which satisfy ( p8[) are known 
as scaling functions. It easily may be verified that if <fi(x) satisfies condition (^), then all basis 
functions 4>(2®y — n) in Vq are also in the space Vq+i, so that indeed Vq C Vq+i- 

Because any multiresolution analysis must be based upon scaling functions, it is useful to have 
a prescription for generating all allowable scaling functions. The two-scale relation (|5^) expresses 
the function (p(x) as a convolution of the discrete sequence c n and the continuous function <j)(2x). 
This convolution takes the familiar form of a product in Fourier space, 

d d x 



<m = /^^(i) (59) 

d d x 



e- ifc '^c n 0(2x 



n) 



J (2vr) d 
= (E ^e^A Hk/2) 
= m (k/2)4>(k/2), 

where we have defined mo, the two-scale symbol for the scaling function (p(x), as 

m o(k) = £ |e- ifc -«. (60) 

n 

By construction, the two-scale symbol is periodic on 2-kCq, the hyper-cubic lattice of lattice constant 
2ir. Note that the first line of ( |59| ) gives the normalization convention which we shall use for Fourier 
transforms throughout this work. 

The recursive nature of Eq. (1591) implies that the Fourier transforms of scaling functions always 



take the form 

0(h) = m (k/2)mo(k/A)...m (k/2 n )...j)(0) 

= (f\m (k/2 n ))m. (61) 



\ra=l 



Conversely, substitution into the final lines of (|59|) verifies that any function constructed from (|6l| ) in 
conjunction with any 2-7rCo-periodic function mo satisfying m-o(O) = 1 (so that the infinite product 



converges) results in a proper scaling function. Eq. (61) therefore places the set of acceptable 
scaling functions as in one-to-one correspondence with 2-7rCo-periodic functions with value unity at 
the origin. The interested reader may wish to consult [Dau92] for a rigorous discussion of technical 
mathematical details related to the infinite product. 

IV-D Two-scale decomposition theorem 

With the allowable scaling functions <j){x) and hence the spaces Vq defined, we now turn to the 
issue of what functions form appropriate bases for the detail spaces Wq+\. 
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Before specifying a basis for Wq+i = Vq+i — Vq, we first determine what functions belong to 
this space by establishing the nature of the excess freedom which exists in Vq+i over that which 
exists in Vq. If F and / are two arbitrary functions from the spaces Vq and Vq+i, respectively, 
then 

F(x) = Yl Fn<f>(2 Q x ~ n) 
nez d 

fix) = E M(2 Q+1 x-n). 

neZ d 

These expressions are both convolutions of form similar to (|5^) and also give products in Fourier 
space: 



F(k) = (E^-W 2e )^(^) 



2Qd 

= m F (k/2 Q )4>(k/2 Q ) 

= m F (k/2 Q )m o (k/2 Q+1 )0(k/2 Q+1 ), (62) 
where we have used the two-scale relation flBl^); and 



/» = fE^^ (fc/2Q+1) ' n V(^ Q+1 : 



2(Q+l)d 

\ lb / 

= m f (k/2 Q+1 )0(k/2 Q+1 ), (63) 

where rap(q) and mf(q) are two new 2-7rCo-periodic functions defined from Fourier series constructed 
with the sequences F n and f n respectively. Thus, F(x) and f(x) are in Vq and Vq+i if and only 
if their Fourier transforms are of the forms in (|62| ) and (|63| ) with 27rCo-periodic functions rap and 

We can use these results to verify quickly that any function F in Vq is also in Vq+i. First note 
that if F is in Vq, then there exists a suitable 2-7rCo-periodic function nip to satisfy (|62|). Using 
this function, then clearly the choice 

m f {q/2)=m F (q)m (q/2), (64) 

will cast F in the form (|63|). All that remains to confirm that this places F in Vq+i is to show 
that the function raj defined as in ( |6"4| ) is indeed always 2-7rCo-periodic, as is easily verified using 
the fact that mo and rap both have 27rCo-periodicity. 

Naturally, it is not possible to describe a general fine-scale function / from Vq+i as a function 
F from the coarser space Vq . To write a function of the form (|63|) in the form ( p2] ) , would require 



rap(q)=ra f {q/2)/ra (q/2), (65) 

to be a 27rCo-periodic function. The function mf(q/2), however, can be any 47rCo-periodic function 
of q. The failure of / to be in Vq thus may be viewed as the inability of the 27rCo-periodic function 
m F{q) to describe the full freedom present in raf(q/2), a general 47rCo-periodic function of q. 

As a hint as to how to capture the this full freedom, we note that the Monkhorst-Pack theory of 
Brillouin zone sampling [MP 76] shows that any 47rCo-periodic function may be expressed as a sum 



of 2 d separate functions, each with a special 27rCo-periodicity. The particular 2-7rCo-periodicities 
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chosen by Monkhorst and Pack are those associated with points of the reciprocal lattice of AttCq 
which fall in the first Brillouin zone of the reciprocal lattice of 2ttCo, namely the set of k-points 
ki =rji/2 where 

m € {0, l} d . (66) 

For instance, in d = 3 dimensions, 

= {(0,0,0), (0,0,1), (0,1,0), (0,1,1), (1,0,0), (1,0,1), (1,1,0), (1,1,1)}. 

Thus, all 47rCo-periodic functions of q may be written as a linear combination of functions Pi(q) 
with the special 27rCo-periodicities described by the condition that 

Pl ( q + Q)=e-^/2) pt ( q) 

for all Q in 2ttCq. For convenience, we will adopt the notation that r/Q = and will index sums 
which omit this zero vector by a, except for two cases noted explicitly in the text in Sec. |V-A . 



From the theory Monkhorst and Pack, we thus conclude that if we take the 27rCo-periodic 
function mF(q), which has the same periodicity as Po(q) above, and augment it with 2 d — 1 new 
functions rriG a (q) with the periodicities of the p a (q), 

m Ga (q) = E Ga,ne" j(ri+W2) - 9 (67) 

neC 

for some set of Fourier coefficients G an , the collection of these functions represents the same 
freedom present in the 47rCo-periodic function mf(q/2). Thus, for any set of fixed 2-7r-periodic 
functions m a (q) introduced to play roles parallel to that of mo(q), we may always find functions 
m G a (q) of the form ( p7|) so that 

2 d -l 

m f(l/ 2 ) = m F (q)m (q/2) + ^ m Ga (q)m a (q/2), (68) 

a=l 



provided the choice of the fixed functions m a (q) is not pathological. Eq. (38), for instance, will 
not hold in general if we choose m a (q) = for all a. For reasons which will become clear in the 
next section, the m a (q) are called the detail or wavelet symbols. 



The Two- Scale Decomposition Theorem [Dau92, Chu92a] specifies precisely the conditions under 



which the choice of m Q (q) is not pathological so that m,f(q/2) indeed may be decomposed in the 
form (|68|). Translating ( p8| ) by each of the 2 d vectors 27r?7j and simplifying the result using the 
periodicity properties of rap and mc , gives the following system of equations: 

m f(q/ 2 ) = m F (q)m (q/2) + ^ m Ga (q)m a {q/2) (69) 

a=l 



2 a -l 



mf(q/2 + irrn) = m F (q)m (q/2 + nrn) + E (- 1 ) Vl ' ri °" m G a (q)m a (q/2 + nrji) 



a=l 



2 a -l 



m MI 2 + 7 ™?2 d -i) = rn F (q)m (q/2 + ir^d^) + ^ (-l)^" 1 rJa rn Ga (q)m a (q/2 + 7rr? 2 d_ 1 ) 



a=l 
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So long as this 2 d x 2 d system is not singular, we may always solve directly for the variable functions 
m f and mc a in terms the function m f we are to reproduce and the 2 d fixed functions mi and thus 
accomplish the two-scale decomposition 



Replacing k = q/2, the final condition for two-scale 



decomposition to obtain is therefore 
m (k) 

m (k + nrfi) 



m 2 d_ 1 (k) 



(-l) n ^ d -^ m m 2d _ l (k + urn) 



m (/c + 7r? ? 2d_ 1 ) ... (-l) V2d - 1 ' V2d - 1 m 2d _ 1 (k + TTrj 2d _ 1 ) 



/ for all k. 



(70) 



We conclude this section by observing the implications of ( f70| ) for functions in the spaces Vq 
and Vq+\. Making the replacement q = k/2® in (p8|) , which we now know holds whenever the 
determinant ( f70| ) is non-zero, multiplying through by (j)(k /2^ +1 ) and using (|62"|), we obtain 



2 a -l 



f(k) = m F {k/2 Q )mQ{k/2 Q+l )4>{k/2 Q+l ) + ^ m Ga {k/2 Q )m a {k/2 Q+1 )(f){k/2 Q+1 ) 

a=l 

= F(k) + d(k), (71) 



which decomposes a general function / G Vq+i into a coarse function i* 1 in Vq and an additional 
function G, which carries the finer scale, detailed information about / and thus lies in the space 
W Q+1 . 

IV-E Wavelets 

To determine the basis functions for Wq+i, note that by its definition (55), the detail space contains 
precisely the functions G(k) from (]7l|). Using the expansions ( |67|) for the functions m^, the Fourier 
transform of (fnl) gives 



nGCo a 



n + 77Q./2 



2Q 



where 



The detail space is therefore 

W Q+1 = span {v> a [2 Q (x 



Mk/2 Q ) = m a {k/2 Q+1 )4>{k/2 Q+1 ) 
n + r/ Q /2 x 



2 Q 



nGC ,«=1...2 d -l 



(72) 
(73) 
(74) 



where the r/ a are defined as the non-zero vectors in ( |6q ) . The new functions ^ a (x) , defined in terms 
of the detail symbols m a from (|68|), when scaled and centered appropriately, make up the basis 
for Wq+i and are thus the detail functions, or wavelets. Note that the extra phase factors in (|67|) 
have ensured naturally that each detail function ip a {x) is associated with a unique point in Dq + i, 
as anticipated in Sec. [IV-A . 

The functional form of the detail functions is revealed by Fourier expanding the 27rCo-periodic 
symbols m a , 



m a (q) = 



2 d 



-iq-n 



(75) 
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and transforming the Q = case of ( |73| ) to real space, 

4>a{x) = ^2 d a ^ n (j)(2x -n), (76) 

neC 

Thus, as Vq © Wq + \ = Vq + \ requires, each detail function from Wq + \ is a linear combination of 
scaling functions from Vq+i. The allowed choices for these linear combinations are those which 
satisfy ( |70| ) with the m a (q) defined through ([75|). 

An alternate route to two-scale decomposition is to note directly from Vq © Wq + \ = Vq + \ that 
the ipa{x) must be formed from linear combinations of the scaling functions from Vq+i. Then, 
to ensure that the two-scale and single scale representations are equivalent, one must show that 
each basis function of Vq+i may be expanded in terms of the basis functions from Vq and TVq+i. 
When Fourier transformed, the system of equations which must be solved to find the corresponding 
expansion coefficients is precisely (^) , leading once again to the determinant condition (|7^) for the 
acceptable choices for the detail functions. 

IV-F Matrix language 

To develop a matrix language useful for the application of multiresolution analysis to physical 
problems, we begin by developing a language for the various descriptions of the functions f(x) in 
Vjv provided by multiresolution analysis. There is the single-scale representation, 

f(x)= J2 (F N ) P H2 N (x-p)), 

where we define Fn as a column vector whose components are the coefficients associated with the 
scaling functions of Vjv. Note that we index the components of Fn by the points p on which the 
basis functions are centered. Alternately, we may represent f(x) using the multiresolution analysis 
which starts from a scale M < N, 

V N = V M © W M +i © ... © W N . 

This gives instead 

TV 

/(*)= E (F NtM ) p <f>(2 M (x-p)) + E E (F mM ) p ^(2 p -\x-p)), (77) 
peC M p=m+i P eDp 

where, regardless of its scale, the expansion coefficient associated with the point p of the multires- 
olution representation spanning scales M through N is indexed as simply (Fn-.m)p- Because the 
expansion coefficient associated with each point depends on the scales present in the multiresolution 
analysis, it is critical to include the subscript N : M on the vector (Fn : m) p in order to specify the 
scales in the multiresolution analysis. Finally, we note that under this convention the single-scale 
coefficients are also given by Fjy = Fn-.n- 

Next, the two-scale relation relation allows us to write both the scaling functions on scale M 
and the detail functions on scale M + 1 as linear combinations of scaling functions on scale M + 1. 
With these replacements, ( p% becomes 

f{x) = E (Fn--m)p E c„0(2 M+1 (x -p) - n)) 
p&C M neCo 
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+ E (F N:M ) pf J2 d a(p'),n^ M+1 (x-p')-n)) 
p'£D M+1 neZ d 
N 

+ E E (F N:M ) p ^2 p -\x-p)) 

P=M+2 p&Dp 

where the a(p') in the coefficients d a ^ p i^ n in the sum for Dm+i are defined so that p' = Po + 
Va(p')/^ M+1 f° r some po in Cm- By collecting terms, this expression rearranges into the expansion 
for f(x) in the (M + l)^ 1 scale representation, 

/(s) ee E {FN:M + i) q ^2 M+ \x-q)) 

qeC M +l 
N 

+ E E (^V : M+l)p^(2 P - 1 (x-p)). (78) 

P=M+2 pGDp 

The linear map which this procedure represents between Fn-.m and Fn-.m+i is most compactly 
represented as a matrix equation, 

Fn-.m+i =Fm+i,mFn-.m- (79) 

The wavelet transform, the recovery of the single-scale representation Fn from the multiresolu- 
tion representation Fn : m, is the result of cascading together these representations of the two-scale 
relation, 

Fn = Fjy : N = (1n,N-i(Zn-1,N-2 ■ ■ ■ (Fm+1,mFn-.m))) (80) 

Note that many authors, particularly in the field of signal processing, refer to the operator In-.m as 
the "inverse" wavelet transform. Because the primary mode of operation in physical calculations 
is to treat the multiscale expansion coefficients as the independent variables, in this work we refer 
to this process of producing a function from its multiscale coefficients as the forward transform, 
which conforms to our usage in Sec. III- 



The transform operator appearing in (|80| ) may be defined as connecting any two scales P > Q, 

Q+i 

Fp-.Q = II X *.*-i> ( 81 ) 
R=P 

where in this expression, as throughout this work, we adopt the convention that non-commutative 
matrix products proceed in order from left to right with the lower index of the product appearing 
leftmost in the product and the upper index appearing rightmost. A direct result of the definition 
(H) is that l M +l:M = F M +i,M- Also > 

lp:Q = lp :R lR:Q for (P > R > Q) . (82) 

Finally, to discuss restricted multiresolution analyses, where functions are selectively removed 
from the basis as described in Sec. Ill- A , we introduce projection matrices which "zero out" coef- 
ficients associated with the points not contained in a given grid G, 
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The product of two such projections is the projection associated with the intersection of the asso- 
ciated grids, and, for disjoint grids, the sum of two such projections is the projection associated 
with the union of the two, 



VGi+Vg-z = T > G 1 \JG 2 +^ > G 1 nG 2 - 



(84) 



Most statements about the sets of points Cq, Dq, Bq defined in Sec. |IV-A| may be derived in terms 
of and translated into identities involving the projections V . Two facts of which we shall make 
frequent use are that the finer grids contain the coarser grids and that the detail points of level P 
are contained only in the grids of level Q > P. In terms of projections, these statements are 



Vc 



min (P, Q) 

Q < P 
V Dp Q>P 



(85) 



respectively. Also, because the linear map represented by Zq+^q replaces scaling functions and 
detail functions on scales Q and Q + 1 with scaling functions on scale Q + 1, it acts independently 
of the coefficients associated with scales beyond Q + 1. Thus, 



T > b q+1 1q+i,q = Iq+i,qPb q+1 = Vb, 



Q+l 



(86) 



V Bases for Mult iresolut ion Analysis 

Even with the restrictions of multiresolution analysis as laid down in the previous section, significant 
freedom remains in the choice of the scaling and detail functions. Motivated by differing aspects of 
the calculations, researchers in electronic structure have employed different multiresolution analyses. 
Ultimately, the optimal approach may be to use different bases for different phases of the calculation. 
The simplicity of working with orthonormal basis sets has led several authors dealing with with 



electronic structure] WC9£ , TW97] to use the wavelets of Daubechies[Dau92|. Section V-A below 
reviews the construction of these bases. 

The primary issue in the choice of bases in problems in the physical sciences, however, is the 
ability of the basis to represent the relevant physical functions, not necessarily the analysis of these 
functions into separate frequency components. While critical in signal analysis, orthonormality is 
not as crucial in electronic structure, as illustrated by the great success of the chemistry community 



with the use of Gaussian bases [FD9C]. Several authors therefore have exploited the freedom of not 
being confined to orthonormal bases to improve the representation of physical problems [ |ACLT95 , 
LAE98] , [Yes97| , |FS94| , |LAE98| , [BN96| , pI98[| . The majority of these a pplicat ions use multiresolution 



analyses based upon the scaling functions of Deslauriers and Dubuc jDD89 l, which were developed 



independently in several fields [|Lem91| , |CS9^ , pon92| , |BS9^ , |SB93l |AU93| , [Uns93| , |ACLT95| , [Lew94 l. 
The great advantage of using these scaling functions is that they function as finite element functions 
and thus provide both good interpolation and transform properties. 

Teter |TeT93|| was the first in the electronic structure community to recognize the advantage of 
adapting the ideas of finite elements to bases with a multilevel structure. A key concept from finite 
element theory is the property of cardinality, the condition that each basis function have value 
zero at every point of the finite element grid except for the one point with which it is associated. 
Unfortunately, it is impossible to maintain cardinality in a multiresolution analysis because smooth 
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basis functions from coarse scales cannot oscillate so as to be zero on the grid points of the finer 
scales. Instead, Teter suggested the construction of bases where cardinality is maintained for the 



points of coarser scales but sacrificed for finer scales. Section V-B gives the precise definition for such 
semicardinal bases, which developed out of this idea, and describes the nature of multiresolution 
analyses which are semicardinal. 

Such semicardinal bases have the remarkable property that expansion coefficients for a function 
may be extracted exactly from the values of the the function on a finite set of points in space. 
This property makes these bases ideal for non-linear problems because non-linear couplings such 
as those in the LDA Lagrangian are computed in real-space on a finite grid and the corresponding 
expansion coefficients must be extracted from the resulting values. This highly desirable exact 
extraction property requires the detail functions to have non-zero integral, and thus prevents such 



bases from being wavelet bases in the technical sense described in [Dau92]. Such semicardinal bases 



were used in [ACLT95, Ari95| , |LAE98|j a nd are described in detail in Sec. V-B below. 



Finally, Goedecker and Ivanov[ pI98|j have suggested that the Poisson equation be solved with 
basis functions with zero integral and higher multipole moments, as formed by linear combinations 
of cardinal scaling functions according to the lifting scheme of S weldens [ |5 we96 ] . Because working 



with functions of zero integral extends the range of the basis functions and disrupts the exact recon- 
struction property which is so useful in treating non-linear interactions, an intriguing possibility for 
future research would be to use semicardinal basis for the non-linear phases of the LDA calculations 
and a lifted basis for the solution of Poisson's equation, should lifted wavelets prove superior in the 
solution of Poisson's equation. (See discussion in Sec. 



V-A Orthonormal bases 

Daubechies was the first to construct orthonormal bases of wavelets with compact support [Dau92] 



We now give a brief review of the discussion in [Dau92] within our present conventions and language 



The scaling functions on the coarsest scale of an orthonormal multiresolution basis are them- 
selves orthonormal. This condition may be re-scaled to the scale Q = 0, where it becomes 

5 nm = [ d d x 4>*{x — n)<f){x — m) 

d d x ( J d d k ^(k)e ik < x ~ n ^)*( J d d k' j>(k , )e ik '< x - m 'i 
(2ir) d / d d k <p*(k)(j)(k)e ik < n - m) for all n,m in C . 



The entire reciprocal space of vectors k may be partitioned into images of the first Brillouin zone 
centered on the points G of the reciprocal lattice 2ttCq. Writing the integral in this way gives, 



5 nm = (2ir) d f d d k 4>*{k + G)4>{k + G)e 

JB.Z. G£2itC 



ik'(n—m) 



where we have used the fact e lG '( n m ) = 1 for the reciprocal lattice vectors G. The fact that this 
is true for all n,m in Co, combined with the Fourier space version of two-scale relation (|59|), gives 



(2ir)- 2d = Y, 4>*(k + G)<f>{k + G) (87) 

Ge27rC 

^ / ,k + G.,,k + G.y{ , k + G ^j., k + G < 

= [ m o(^—m^—)j [mo{— — M— =— , 
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Finally, defining q = k/2 and re-expressing the sum over G/2 as a sum over the reciprocal lattice 
and the decoration points nrji with 77, defined in (|6lf) , gives 



(27T) 



-2d 



E E m o(</ + (G" + tt Vi ))0* (q + (C + ir Vl ))m (q + (G' + vr^))^ + (C + tt^)) 

G'e27rC i 

E™o(<7 + ^W<Z + ^) E ^*(« + i G ' + ^))0(« + ( G ' + 
j G'e27rC 

E m o(<? + Trrji)m (q + vr77 i )(2vr)~ 2d , 



where we have used the periodicity of rriQ and the first line of (|37|). We thus conclude that or- 
thonormality among the scaling functions places the following condition on the two-scale symbols, 

2 d -l 



E m o(<7 + ^r]i)m (q + nrji) = 1. 



i=0 



In Sec. V-C , we show how to construct functions M^{q) which obey the constraint 



2 a -l 



E M (g + 7nfc) = 1. 



i=0 



Comparing with 



we see that to construct orthonormal scaling functions, one first finds such 



an Mo(q) and then takes the "square root" to find the appropriate two scale symbol mo(q). This 



procedure is described in Dau92 , Sec. 6.1], which also tabulates the resulting two-scale coefficients. 

Next, to construct orthonormal detail functions, one must assure that the wavelet tp a for each 
decoration point a = 1 ... 2* — 1 be orthogonal to the scaling functions, 







d d x (j)*(x — n)ip a (x — m — r? a /2) 



ik-(m—n)—i(k+G)-ri a /2 



= (2vr) d / d d k Y$*(k + G)Tp a {k + Gy 
Jb.z. q 

Through an analogous analysis to that above, we will finally arrive at 

= E m o(9 + m) (e~ mri ^ a m a {q + nrji)) ■ 

i 

One must also ensure orthonormality among the i/j a (x), 

Saturn = J d d x $*Jx - n- r] a /2)ip /3 (x - m - 77/3/2), 



(89) 



which leads to 



2 fl -l 



i=0 



(90) 



Finally, we note that three orthonormality conditions (88-90) may be combined into one single 
condition by simply taking a and f3 in (^0|) to vary over the full range . . . 2 d — 1. 

In one dimension, Eq. (|90| ) represents three independent conditions: orthonormality among the 
scaling functions, orthonormality among the wavelets ip\{x — n — 1/2), and orthogonality between 
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these two sets of functions. Given a two-scale symbol mo satisfying the orthonormality condition, 
one common choice for the wavelet symbol m\ which satisfies the remaining conditions in (^0|) is 

mi{k) = mo(fc + vr), (91) 

as may be confirmed directly by substitution. 

For d > 1 dimensions, a straight-forward way to satisfy the conditions ( |90| ) is to use product 
functions formed from the symbols mo and m\ of an orthonormal, one-dimensional multiresolution 
analysis, 

d 

m a (q) = ]J m (^) e (^), 

e=l 

where e varies over the Euclidean dimensions of space, q e and (r] a ) e are the corresponding compo- 
nents of the respective vectors, and, again, we let a vary over the entire range including a = 0. 
With this choice, the condition ( |90| ) factors and gives 

2 d -l 
i=0 



= IL [Ej=o (^^""H^te + nr,))* (e-™^m (m)e {q + irq)) 
= n e *(j7a) e ,(^)e = S a,/3, 

as is required. This, coupled with the choice (|9~l]), was used to construct the multidimensional basis 
functions used in the electronic structure calculations of [ WC96 , TW97]. The particular choice of 
two-scale symbols tuq in these works were the D6 and D8 wavelets defined in [ Dau92j 1 , respectively. 



V-B Semicardinal bases 

It is highly desirable in the calculation of electronic structure to be able to determine the expansion 
coefficients of functions from knowledge of their values on a grid of finite resolution. If we insist 
that at every level of resolution Q, the expansion coefficients of a function depend only upon its 
values on the grid Cq, then, as we will now see, it is always possible to find a multiresolution basis 
in which every basis function has the value zero on all points of the grid of its corresponding scale, 
except for one, where its value is unity. We refer to such a basis as a semicardinal basis. 

The following definitions are useful in discussing such bases. The concept of semicardinality 
may be defined independently of multiresolution analysis, so we couch these definitions in general 
terms: 

1. A hierarchical basis is a set of functions, each associated with exactly one point g from a 
hierarchy of nested grids Gm C Gm+i C . . .. 

2. For such a set of grids, Hq + \ = Gq + \ — Gq. 

3. In a hierarchical basis, the scale of a point g and it associated basis function is the least value 
of Q for which g is in Gq. 

4. A function f(x) is cardinal on a given grid G if f(g) = 5 g h for all points g and one point h, 
both in G. 
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5. A semicardinal basis is a hierarchical basis in which each basis function is cardinal on the 
grid of its own scale. 



Note from the second and third definitions that when Q = M, Gq is the set of all points from scale 



Q and that when Q > M, Hq is the scale of all points from scale Q. In the language of Sees. III-C 



and |IV-F , semicardinality is the condition that the matrix I, whose columns contains the values of 



each basis function on the finest grid, obeys two simple algebraic identities, 

Tg m 1V Gm = V Gm (92) 
Vg q 1Vh q = V Hq for Q>M. 

Note that a multiresolution analysis is a hierarchical basis on the grids Cq , with the Dq playing 
the role of the Hq. As an example, consider the multiresolution basis consisting of the basis 
functions from the top two panels of Figure ||. The basis functions for Vo are associated only with 
points of Co but not C\ and thus are of scale N = 0. These functions are cardinal on the grid 
Co and thus satisfy their condition for semicardinality. The basis functions for W\ are associated 
with the points in D\ = C\ — Co, and thus are of scale N = 1. But, they are not cardinal on the 
grid Ci, and thus this two-scale basis fails to be semicardinal. For the basis to be semicardinal, we 
would have to replace the functions centered on the odd half-integer points with functions which 
are cardinal on C\. 



V-B.l Exact extraction 

In a hierarchical basis on the scales N : M (by our convention N > M), the property of exact 
extraction means that for all scales Q, the expansion coefficients on the scales Q : M for a function 
f(x) may be determined from knowledge of the values /(Gq) of the function on the grid Gq only. 
This means, in particular that the map from the expansion coefficients on scales Q : M to the 
values of the function /(Gq) is invertible. We now show that this property of obvious practical 
utility implies that such a basis may always be chosen to be semicardinal. 

Consider functions expanded from the coarsest scale M up to a very fine scale R » Q. By 
our assumption of the existence of an invertible linear map, the expansion coefficients on the scales 
Q : M of any function with values /(Gq) = must be zero. Thus, the space of functions satisfying 
/(Gq) = must be a subspace of span of the basis functions associated with the remaining points, 
Gr — Gq. 

Now, the space of functions /(Gq) = has exactly one degree of freedom for each point of 
Gr — Gq and thus has the same dimension as the span of the basis functions associated with the 
same set of points. Having already established that the former space is contained in the latter, we 
now see from their dimensionality that the two spaces are in fact the same. 

Specifically, this means that all functions f(x) in the span of the basis functions associated with 
Hq + \ must all have zero value on the next coarser grid, /(Gq) = 0. Hence, the invertible linear 
map which we have supposed to exit between the coefficients of scales Q + 1 : M and the values 
/(Gq+i) induces an invertible map from the span of the basis functions associated with the points 
Hq + \ to the sequences of values /(Gq + i) for which /(Gq) = 0. A basis for this latter sequence 
space is the set of sequences which are zero on all points of Gq+i but for one point in Hq + \ where 
they take the value unity. Applying to each of these sequences the inverse of the aforementioned 
induced map produces a basis for the span of the functions of scale Q + 1 satisfying the conditions 
of semicardinality. Proceeding in this manner for all scales N < Q < M produces a semicardinal 
basis. 
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Figure 10: Sparsity pattern of Hermitian operator in a semicardinal basis on Cq,Ci,... in one 
dimension. (Basis functions grouped into blocks according to scale from coarsest to finest.) 



V-B.2 Algorithms 

Although semicardinal bases have the desirable property of exact extraction, without the additional 
structure of multiresolution analysis, it is difficult to develop efficient methods for performing the 
operations needed in physical calculations. 

In particular, without the structure of multiresolution analysis, the operations C, 0,I,1> must 
be applied directly as multiplication by the corresponding matrices. When working with the com- 



plete grids Co,C\, . . . as defined in (p4|), these matrices show a sparse, fractal- like pattern. Figure 10 
shows the appearance of this pattern in one dimension for L and O when the basis functions are 
ordered first according to scale from coarsest to finest and then within each scale by location in 
space. A similar pattern results for the matrices X and except that, because of semicardinality, 
these matrices are either zero above or below the diagonal, respectively. 



After the process of restriction, as defined in Sec. III-A and illustrated in Figure [|, little of 
this sparsity remains. In calculating electronic structure, the majority of functions surviving the 
restriction will overlap with one of the nuclei. The matrix elements among these surviving functions 
thus consist of dense blocks connecting all of the functions associated with each atomic nucleus. 
Multiplication by matrix blocks of such size requires thousands of operations per basis function 
and is therefore relatively inefficient. 

The inverse transform J has the potential to require the solution of a general system of linear 
system of equations and thus become even more problematic. Semicardinality, however, provides 
just enough structure to allow the inverse transform to be performed in a direct procedure requiring 
the same work as does X. To see this, let / = f(Gjsr) be a column vector of the samples of a 
function on the finest grid and F = J f be the expansion coefficients we seek. Decomposing F 
into contributions on different scales and applying the identity Vg m TVh q = 'PgmCPgqZ'Phq) = 
Vg m 7~ > Hq = 0, which follows from the algebraic definition of semicardinality (|92| ) and the facts that 
Gq C Gm and Gm H Hq for Q > M, gives the following result for the values of the function on 
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the coarsest scale, 



V GM f = V Gm IF (93) 

= v Gm i(v Gm +v Hm+1 + ...)f 
= ?g m f. 

Thus, the expansion coefficients on the coarsest scale are just the values of the function on the 
associated points. 

Now, proceeding iteratively, suppose the coefficients of F are known up to scale Q and consider 
the values of / on the next scale. Again using the conditions (|92|), we have 



Vh q+ J = V Hq+1 1{V Gq +V Hq+1 +V Hq+2 + ...)F 
= (V Hq+1 1V Gq +V Hq+1 )F, 

which gives the expansion coefficients on the next scale explicitly as 

V Hq+1 F = Vh q+ J-V Hq+1 FV Gq F (94) 



Combining (^,94) gives the following recursive procedure for computing F = J~f, 

F M = /; (95) 
F Q+ i = F q -V Hq+1 FV Gq F q ; 
F = F N . 

Note that the total action of all of the operators T , Hq +1 ZT , Gq is just to connect once each of the 
columns of I with all points in the Hq + \. The implementation of (^) is thus precisely as costly 
as the implementation of X. 

Finally, note that the algorithm ( |95"|) may be written as the matrix product, 



M 

J= II (i--Ph q+1 1V Gq ), (96) 

Q=N-1 

from which the procedure to compute is easily determined. 
V-B.3 Semicardinal multiresolution analysis 

The key feature of semicardinality, that the expansion coefficients of scale Q for a function f(x) may 
be extracted from the sample values /(Cq), is incompatible with the use of standard orthogonal 



or bi-orthogonal multiresolution analyses. In bi-orthogonal multiresolution analyses | Dau92 1 , one 
considers both the basis of scaling and detail functions and the dual to this basis. The dual to any 
basis is defined as those functions against which any function f(x) may be integrated to determine 
its expansion coefficients in the original basis. Orthogonal multiresolution analyses are thus a 
special case of bi-orthogonal bases where the basis is its own dual. For semicardinal bases, the 
exact extraction property implies that the dual of the functions of scale Q are linear combinations 
of Dirac 5-functions centered on the points of Cq. Thus, a semicardinal basis of smooth functions 
can never be orthogonal, and, because Dirac 5-functions are not square-integrable, such a basis 
does not technically fit into the bi-orthogonal wavelet framework. 

Nonetheless, it is possible to construct multiresolution analyses which are semicardinal. Semi- 
cardinality, in fact, nearly completely determines the allowable form for a multiresolution analysis. 
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To form a semicardinal multiresolution analysis it is necessary and sufficient that the scaling func- 
tions be cardinal and that each of the 2 d — 1 detail functions for scale Q simply be the scaling 
functions of scale Q + 1 associated with the corresponding detail points. While this latter property 
may seem unusual at first, having only one type of function in the multiresolution analysis proves 
to be of considerable convenience |LAE98| . 

To show that semicardinal multiresolution analyses must have this form, consider first the 
requirements which semicardinality places on the scaling functions. The first semicardinality con- 
dition of ( |92| ) states that the functions spanning the coarsest scale M must appear cardinal on the 
coarsest grid, which is Cm in the case of a multiresolution analysis. Dilating this condition to the 
scale Q = 0, gives the condition that the scaling functions must be cardinal on the integer grid, 



<j)(n — m) 



for n, m in Cq = Z a 



(97) 



Section V-C.l(b) discusses the implications of this for functions satisfying the two-scale relation. 

Next, we turn to the detail functions. The second semicardinality condition of (]9^ ) states 
that for all q E Dq and p € Cq, the detail function associated with point q must have value 
tp a (2 ( ^~ 1 (p — q)) = 5 pq . Enforcing this condition for all p and dilating the condition to scale Q = 0, 
implies that for all a, ip a (n/2) 
of (f7§) that 



<Ln for n G Z . From this we find that for the detail coefficients 



kez d kez d 
= ip a (n/2) = 5 n0 . 



(98) 



Therefore, as stated above, the detail functions tp a (x) = J2k&z d d a ,k<ft(2x — k) = 4>{2x) are just 
the scaling functions from the next finer scale. For a general scale Q this means that the detail 
functions for the points Dq + \ are just the scaling functions on the corresponding points in the 
single-scale representation of Vq + \. 

Finally, we now verify that such a basis satisfies the two-scale decomposition theorem. Eq. 
implies that m a (q) = l/2 d for all a and q. Thus the determinant from (l70|) becomes 



det 



/_|_\2 d -l 

y 2 d> 



m (q) 



m Q (q + nrn) 



m (q + vr77 2 d 



-1/ 



(_1 ^Va-Vi 



1 



(_l)»?ad-i-V-i 



(99) 



where, we have used the fact, derived in the appendix, that each cofactor in the determinant 
expansion along the first column has the value ±(2 rf ) 2 _1 , with a fixed sign for each cofactor. For 
cardinal scaling functions, we have from ( \10l 



in Sec. V-CT , that the sum over the mo appearing 
in (^9|) is unity. Thus, as required for two-scale decomposition, the determinant is never zero. The 
fact that we find a simple constant for the determinant is a direct reflection of the compact nature 
of the dual basis. 
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V-C Interpolating, cardinal scaling functions 

Because the previous section has already determined the form which the detail functions of a 
semicardinal multiresolution analysis must take, it remains now only to specify the scaling func- 
tions of such multiresolution analyses. Accurate calculations in physical systems require that the 
reconstructions of physical functions from sample values interpolate well the behavior of those 
functions between the sampling points. Additionally, in order to limit the processing required for 
each expansion coefficient, in practice the basis functions should have the minimal spatial extent 
possible. These conditions represent two additional constraints beyond those we have already im- 
posed, cardinality for exact recovery of expansion coefficients and the two-scale relation to sustain 
multiresolution analysis. As we shall show in Sec. [V-C.l , the four constraints of (1) the two-scale 



relation, (2) cardinality, (3) minimal support, and (4) interpolation almost uniquely specify the 
scaling functions. In the interest of brevity, we shall refer to such interpolating, cardinal scaling 
functions as "interpolets." 

Once the two-scale coefficients of the interpolets are determined, the only other information 
needed to perform physical calculations are the matrix elements of integral-differential operators 
between scaling functions and the values of the scaling functions in real-space. Sections V-C.2| and 



V-C. 3 describe the determination of these quantities, and Sec. V-C. 4 gives examples of specific 



functions, with tables displaying all information needed in physical calculations. 
V-C.l Construction 

We now discuss the mathematical implications of the four constraints outlined above. Because we 
now consider the very specific class of interpolating cardinal scaling functions, we shall denote these 
functions as I(x), in place of the generic 4>(x). 



V-C. 1(a) Two-scale relation As discussed in Sec. [IV- C , we may consider the entire allowable 



class of scaling functions by considering only functions constructed in Fourier space via 

^)=|fi^'m(o) ( 10 °) 

with an arbitrary 2-7rCo-periodic function mo(q) satisfying mo(0) = 1. 

Each of the remaining conditions thus becomes a condition on the acceptable two-scale symbols 
mo(q). Once mo(q) is known, the two-scale coefficients c n 

l{x) = Y,c n 1{2x-n), (101) 



used to form the operators Zp+i p from Sec. [IV-F| , may be recovered from (J60j 



V-C. 1(b) Cardinality Cardinality is the condition that the scaling functions X{x) obey 

J( n ) = S n>0 , for all n £ C . (102) 

To transform this into a constraint on the c n , and thus on the mo(q), we dilate the two-scale 
relation by a factor of two, 

X(x/2) = 5> n J(x-n), (103) 
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and evaluate the result on the points x = m in Co, 

J(m/2) = ^2 c n1{m ~n) =Y1 C n^mn = Cm- (104) 
n n 

, cardinality determines all "even" terms in Cn, 

C2n = Ko, for n in C . 

To convert this into a condition on the two-scale symbol, we remove all "odd" frequency com- 
ponents of mo (q) leaving untouched the "even" frequency components by periodizing the sum ( |60|) 
on the lattice ttCq- Cardinality thus implies that the two-scale symbol obey 

2 f^mo(k + 7r Vi )= J2 2^e- lfc - = l, (105) 

i=0 n62C*o 

which we used in ( |99|) to confirm that semicardinal multiresolution analyses indeed satisfy the 
two-scale decomposition theorem. 

V-C.l(c) Minimal support From ( p.04|) we see that cardinal scaling functions extend at least 
as far as the range of non-zero elements in the sequence c n . We thus limit our discussion to 
sequences with the shortest possible length. Correspondingly, the mo (q) are to be constructed as 
the finite trigonometric polynomial of the lowest order which satisfies all other conditions on the 
scaling functions. 



Comparing with (102 



V-C.l(d) Interpolation By ensuring the equivalence of the multiscale and single-scale repre- 
sentations, multiresolution analysis simplifies the consideration of how well a multiresolution basis 
interpolates physical functions to the consideration of interpolation for the single-scale represen- 
tation on the finest scale of the analysis. On the finest scale, cardinality of the scaling functions 
gives 

f N (x) = f(p)l(2 N (x-p)) (106) 

peCjv 

as the function f(x) in Vn which matches exactly the sample values /(Cn)- To ensure that this 
estimate interpolate fix) to order I between the sample points, we insist that it reproduce exactly 
all polynomials up to degree I. By linearity, to do this, we need only impose the reconstruction of 
all multinomials Ili^! 1 where — I an d < Zj. Rescaling this condition to the scale Q = 0, it 
becomes 

U x i = ECTN^-n). (107) 

i neCo i 

To determine the constraints which this condition places on mo{q), we first consider the con- 
straints which it places on T{k). Condition (107) is a convolution leading to the familiar product 
form in Fourier space 

-ik-x 

d {li} S^(k) = l(k).(27T) d d {h} S^(k -27m), (108) 
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Figure 11: Fourier transform of an interpolating, cardinal scaling function. 

where 5^ d \q) is the d-dimensional Dirac <5-function and dn.y denotes the mixed partial derivative 
rii^x j an d we have used the identity J2n e ~ tk n = 2~2n S^(k — 2irn). Eliminating the sin- 

gularities at non-zero points in 2ttCq from the right-hand side of (|10§| ) requires a high degree of 
regularity in 1{k) at these points. Integrating the derivatives of the <5-functions by parts reveals 
that ( |108| ) implies the conditions 

(2vr) d 5 {iJ J(2vrn) = \J[ 5 hfi ^ 5 nfi for h < I- (109) 

Figure [ll| illustrates these conditions on I(k) for one dimension. One important result is that 
interpolation sets the normalization of the of the scaling functions. The / = condition at n = 
is just T(0) = l/(27r) d , or equivalently, 

J d d x T{x) = 1. (110) 
The other conditions at n = are that the higher order integral moments up to order I be zero, 

J d d x {\{x\ l )l{x) = for Y,ik < I- 

i 

The conditions at the non-zero points in 2ttCq are most simply expressed as the condition that 
T{k) have an Z^-order ze ro at these points. 
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Figure 12: Two-scale symbol tjiq for an interpolating, cardinal scaling function. 



The n = conditions of ( |109| ) constrain the behavior of mo(q) near q = 0. With a finite number 
of non-zero c n from the condition of minimal support, we are assured that mo(q) is analytic. 
Recalling that mo(0) = 1, we have that near q = 0, 

mo(g) = l + 0(/), (111) 

for some leading order j3. Eq. (100) then ensures that I(k) has a similar analytic structure near 
k = 0, 



j(*o = j(o) n (i + o ((fe/2 



n\/3 



n=l 



J(0) (l + O (*/)) . (112) 



Thus, we satisfy all n = conditions from Q10S| ), including the zeroth-order normalization condition, 
so long as we take 

-I oo 

i{k)= (2^n™o(W, (113) 

for some mo satisfying 

% } m (0)=n^,o forallE^i<Z- (114) 

i 

The final conditions concern T(k) at the non-zero points of IixCq. From the form of the product 
(|113| ), we see that if i(k) is to have an order zero at G, then the orders of the zeros of mo(q) 
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among the points q = G/2, G/4, . . . must sum to give at least this order. Because all non-zero points 
in ttCq = 27rCo/2 may be expressed in terms of the product of some non-negative power of two with 
some vector in ttDi, if mo(q) has an fi 1 order zero at each of the points vrl?i, then I(k) will have 
the appropriate analytic structure at the non-zero points of 2ttCq. Noting the 27rCo-periodicity of 
the mo(q), the final conditions imposed by interpolation on the two scale symbol are thus 



dy^m^irrja) = for all ^.Ji <l and 1 < a < 2 d — 1, 



(115) 



with the r\ a defined as in (pq). Figure 12 illustrates these conditions for one dimension. 



An interesting alternative to ( 115 ) is to note that one could ensure the proper behavior in 1{k) 
by instead placing fi 1 order zeros in mo(q) at the points of 7rL>2- But, this and other such choices 
demand higher degrees of oscillation in mo(q) and would to require longer non-zero sequences in 
c n . Also, condition (115) has the distinct advantage that when combined with condition (|11 



the 



two conditions automatically imply cardinality, condition ( |105| ) 

To minimize the support of the resulting functions, we thus consider scaling functions con- 
structed according to (|113| ) with irio(q) satisfying the conditions ( |114 , 115| ) and which therefore 
satisfy all of our constraints. For a given degree of interpolation I, the conditions ( |114 , 115| ) are 

equivalent to a system of simple linear equations for the two-scale coefficients c n . The most com- 

th 

pact set of coefficients satisfying these conditions gives the most compact I -order interpolating 



cardinal scaling function, or "interpolet" . Section [V-C.4| details these equations and their solution 
in several useful cases. 



V-C.2 Overlaps 

We now present a method for determining matrix elements among the scaling functions directly 
from the two-scale coefficients c n . As many of these results are general, we revert temporarily to 
the notation <p(x) for the scaling functions. Once we begin to use specific properties of interpolating 
cardinal functions, we shall return to the notation I{x). The approach below has been developed 
by several authors in one dimension and is described in [ SN96| |. We here give the appropriate 
generalizations for multiple dimensions. 

The specific information needed to apply the overlap matrices O and L when using the proce- 



dures of Sec. VI is the set of matrix elements among the functions of the finest scale, 

M pq = [ d d x^(2 N (x-p))Mc/ )m (2 N (x - q)) for p,q in C N , (116) 



with M being the the identity operator 1 (for O) and the Laplacian V 2 (for L). To compute these 
elements we exploit the fact that both 1 and V 2 are homogeneous operators. 

The homogeneous operators of order h are those which may be written in Fourier space as a 
linear combination of multinomials of order h, 



M = MW(fc)= ]T m {hj} Y[k^ (117) 



T, i= i h J= h 



j 



for some set of coefficients m^.y. In real space, these operators are just 



M = ( - ) , (118) 



51 



where i = Taken together, the operators of order h thus form a linear space of dimension 

equal to the corresponding number of multinomials, {d — l + h)\/[{d— l)!/i!]. The multinomial 
coefficients m{hj} f° r an y such operator may be extracted simply by acting with the operator on 
the multinomial functions of order h, 

mM = ~iw ■ (119) 

Finally, we note that the identity operator M(k) = 1 belongs to the one-dimensional space of 
operators of order h = and that the Laplacian operator M(k) = — Ylj=i k] comes from the space 
of operators of order h = 2, which in d = 3 dimensions has dimension six. 

For any such homogeneous operator, the required matrix elements ( |116| ) for any scale may be 
related to a single universal set of matrix elements, 

M pq = I d d x 4>*(2 N (x 



J d d x 0*(2 N (x-p))2 Nh [M^ {^\ <j>{u)\ 



u=2 N (x-q) 

= 2 N ( h ~ d \M) 2N{p _ q) , (120) 

where we have changed integration variables to give an expression in terms of overlaps on scale 
Q = only, and where 

{M)neC = / d d x <j)*(x - n)M4>{x). (121) 



Now, to determine M in terms of the two-scale coefficients, we generate a recursion by using 
the two-scale relation to write <j)(x) as a sum of scaling functions from scale Q = 1 and then using 
( |12C| ) to re-express the resulting set of overlaps back in terms of M, 

(M) peCo = I ' d d x4>*(x-p)M<j>{x) 



( d d x ( £ c* m 4>*(2(x - p) - m) I M I c nH2x ~ n) ] 



— \ " „* „ nh-d jTf 

— / y ^rn^n^ lvl 2p+m—ni 

m,n£Co 

which can be rearranged as 

2- h (M) peCo = E \ 2 ~ d E C U P+ nCn) {M) q . 

q&C \ n€Co / 

Thus, for each homogeneous operator of order h for which the vector of matrix elements M is 
well-defined, M is an eigenvector of eigenvalue A = 2~ h of the two-scale matrix 

M pq = 2- d E c*^ 2p+n c n . (122) 

neC 

We therefore expect the eigenspectrum of M to consist of clusters of degenerate eigenvalues of value 
2~ h , with degeneracies equal to the dimensionality of the space of operators of order h and with 
eigenvectors containing the overlaps for /i^ n order operators. 
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Determining the overlaps for a specific operator requires additional information in order to 
select the appropriate vector from the corresponding degenerate subspace of M. In the case of 
interpolating cardinal scaling functions this information comes from the integral normalization of 
the scaling functions (110) and the interpolation condition (|107|). To ensure that a given eigenvec- 



tor represents the correct operator, we verify the values of the multinomial coefficients using the 
extraction formula ( |119| ), 

m {hj} = Jd d xl*(x)m {hj} (123) 
d d xT(x) y 3 ' 



J d d x l*(x) ■ f £ (j[nf j - ") j 



hi 



Thus, to determine the vector of overlaps M for an operator of order h with coefficients m{hj} 
one simply forms the linear combination of eigenvectors of ( |122| ) with eigenvalue 2~ h that satisfies 



the condition s (fL23| ). Because the space of such vectors has one dimension for each multinomial 
condition in ( |123|) , this completely determines the vector M. 

V-C.3 Real-space values 

Two approaches are available to compute the values of the interpolets in real space from their 
two-scale coefficients c n . The traditional approach, appropriate for any type of scaling function, is 



to note that the two-scale relation (58), when evaluated at the points x = m/2 + in Cp+i, gives 



<f>(m/2 p+1 ) = <^((m - 2 p n)/2 p ), (124) 

neCo 

a direct formula for the values </>(Cp+i) in terms of the values <j>(Cp). Using this, one may proceed 
iteratively to compute the 4>(Cq) at any desired level of resolution. To determine the values (j)(Co) 
needed to initialize the process, one may generate a self-consistency relation by evaluating the 
two-scale relation (|58|) on the points x = m of Co, 

(f)(m) = c n <fi(2m — n) 

neCo 

= C 2 m -q<P{q)- 

The sequence <fi(Co) thus may be identified as the eigenvector with eigenvalue unity of the matrix 
C m q = C2m-q- For cardinal scaling functions, of course, the cardinality condition ( |102| ) prescribes 
explicitly that the 1{Cq) be a discrete Kronecker 5-function centered on the origin. 

For cardinal scaling functions, a more convenient approach exists which does not make large 



strides through the data-set as does ( 124 ). Because Vq C Vp, I{x) may be written exactly in terms 



of scaling functions on scale P. Moreover, from the interpolation property ( |106|) , the proper linear 
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combination for doing this is just 

l{x)= J2 I(p)l(2 P (x-p)). (125) 
Letting p = m/2 p , evaluating this relation on the points x = n/2 p+l where n 6 Co, and using 



(104) to relate the values 1{C\) to the c n gives 

l(n/2 p+1 ) = J2 c n . 2m X{ml2 p ), (126) 

meCo 

which again may be used iteratively to generate any Z(Cq) from the known values Z(Cq). A variant 



of this approach exists for orthonormal functions and is described in |Dau92| , Ch. 6.5]. 



V-C.4 Examples 

The simplest interpolating cardinal scaling function we shall consider is the first-order interpo- 
let in one dimension. For I = 1 and d = 1, the conditions ( |114 , 115| ) on the two-scale symbol 



are just mo(0) = 1, m,o(vr) = and m' (0) = m' (ir) = 0. The general appearance of such a 
function (Figure |i~2| ) leads quickly to the ansatz, which may be easily verified, that the choice 
mo(q) = cos 2 (2<7) = (1 + cos(q))/2 satisfies these conditions. From this two-scale symbol, we have 
immediately the two scale coefficients are cq = 1, c±\ = 1/2, and c n = for \n\ > 1, as listed in 
Table @. 

Using ( |126| ) to compute the real-space values of T{x), we see that for this interpolet the value 
on each detail point of -Dp+i is given simply by the average of the values on the neighboring points 
of Cp. Using the notation {p} to indicate a complete shell of points related by cubic symmetry, 
so that simply {p} = ±p in this one-dimensional case, Table [l] gives the results of this process. 
As evident from the table, this first-order interpolating function is the piece-wise linear "triangle" 
scaling function sketched in Figure ||. Finally, Table |l] also lists the vectors of matrix elements 
(AfW) n computed according to the prescription ( |121| ) of the previous section. For the values with 

n < 0, note that Mz^ = {—l) h M^ in this and all cases below. 



Figure ( |l3|) illustrates the estimates /q(x) for the function f(x) = sin(27rx) afforded by these 
first-order functions for scales Q = 2 and Q = 3. The bases for Vq provide simple linear inter- 
polation between the sample points /(Cq), and as required by the interpolation condition (107), 



clearly will reproduce exactly the constant and linear functions. As a quantitative illustration of 



the linear nature of this interpolation, Figure 14 shows on a logarithmic plot the root mean-square 
error in reproducing sin(2-7rx) as a function of the sampling rate 2**, as Q varies from one through 
eight. As expected for linear interpolation, to leading order, the error is second-order and the data 
fall along a line of slope —2 in the plot. 

To construct higher-order interpolets, we note that in one dimension Eq. ( |105| ) ensures that the 
conditions on mo(q) at q = tt will be satisfied automatically whenever the conditions ( |l 14| ) at q = 
are satisfied. In terms of the two-scale coefficients, these latter conditions are 

- n ° c n = <5a,o for < a < I. (127) 

For second-order scaling functions this gives three independent linear equations, 

C-l + 1 + Cl + C3 = 2 
— c_i + ci + 3c 3 = 
c_i + ci + 3 2 c 3 = 0. 
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n 





{1} 




1 


1 

2 



X 







{*} 




{1} 


X(C ) 


1 











T(D 1 ) 






1 

2 






1(D 2 ) 




3 
4 




1 

4 




Z(C 2 ) 


1 


3 

4 


1 

2 


1 

4 






n 





1 




2 
3 


l 
fi 







l 

2 




-2 


1 



Table 1: First-order interpolating cardinal scaling function in one dimension: non-zero two-scale 
coefficients (c n ), function values (Z), and overlaps of d h /dx h (Mi h \ as denned in (|121|) ), 




Figure 13: Interpolation with first-order interpolating scaling functions: sin(2-7rx) (solid curve), 
representation in V 2 and V3 (dashed and dot-dashed curves, respectively). 
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Figure 14: Root mean square error in reconstructing the sine function as a function of the number 
of samples per period: first- (solid line), second- (dash-dotted line) and third- (dashed line) order 
interpolating scaling functions, exhibiting exponents of -2, -3 and -4, respectively. 
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-1 
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3 
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1 


3 
4 


i 

8 
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3 
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2 


b 

2 


3 


J(d) 





3 

8 


1 


3 
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1 

8 












n 





1 


2 


3 


MW 


247 
295 


517 
4720 


17 

590 


3 

4720 







11 
16 


1 

10 


1 

240 


M< 2 > 


30 
11 


397 
264 


5 

33 


1 

88 



Table 2: Second-order interpolet in one dimension: non-zero two-scale coefficients (c n ), function 
values (J), and overlaps of d h /dx h (Mn h \ as defined in ( |121[) ). 



n 
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{3} 




1 


9 
16 


i 

16 
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{§} 


{2} 


m 


{3} 


T(C X ) 


1 


9 

16 





1 

16 












n 





1 


2 


3 


4 


5 




56264 


19253 


2827 


6283 


16 


l 


70245 


140490 


70245 


2247840 


210735 


6743520 


M« 





3659 
5280 


731 
6930 


481 
73920 


4 

10395 


1 

665280 


M< 2 > 


20 
9 


9 
8 





1 

72 









Table 3: Third order interpolating cardinal scaling function in one dimension: non-zero two-scale 
coefficients (c n ), function values (Z), and overlaps of d h /dx h (Mn h \ as defined in ( |121| )), 



Table ^ lists the two-scale coefficients which solve this system and the values of the corresponding 
scaling function on the points of C\. Figure ^ shows the appearance of the second-order interpolet 
in real-space. The asymmetry of the function is a direct result of the lack of symmetry in the 
sequence c n . Figure 14 shows the reconstruction of the sine function for these interpolets also. The 
reconstruction is correct to second order and the leading-order error is third-order, as expected. 

The final function we consider in one dimension is the third-order interpolet, which was used in 
the calculations in [ ACLT95 , |Ari95 , LAE98 | to produce the results shown in Figures [j], Hi H and[B|. 
For odd orders of interpolation, the odd-order conditions from (127) are always equivalent to the 
condition that the interpolet be symmetric, c n = c_ n . The remaining conditions for the third-order 
case are just 



ci + c 3 
ci + 9c 3 



1 

2 
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Figure 15: Second-order interpolet in one dimension. 




n 


000 


{001} 


{011} 


{111} 


Cn 
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Table 4: Two-scale coefficients for first-order interpolating cardinal scaling function in d = 3 
dimensions (product form). 



n 


{000} 


{001} 


{011} 


{111} 


{003} 


{013} 


{113} 


{033} 


{133} 


{333} 


c n (3d) 


1 


y 




11 


l 


l 


l 











C n (Pd) 


1 


— tf— 

16 


—to— 

256 


4096 


— ¥- 

16 


256 


— H- 

4096 


l 

256 


9 

4096 


l 

4096 



Table 5: Two-scale coefficients for third-order interpolating cardinal scaling functions in d = 3 
dimensions: full d = 3-dimensional construction (3d), product form (pd). 



Table [3| gives the resulting two-scale coefficients, function values on Ci, and vectors of matrix 
elements for homogeneous integral-differential operators. Figure [l6| shows the function in real 
space, and Figure 14 shows that the leading error in reconstructing the sine function has the 
expected fourth-order behavior. In one dimension, for each odd order, the minimally supported 
interpolets are just the scaling functions of Deslauriers Dubuc |DD89|. 

In higher dimensions, the simplest way to construct interpolets is to take the two-scale sym- 
bol to be an outer product mo(q) = Ilf=i m o fai) °^ ^ ne two-scale symbol ml' (q) for the one- 
dimensional interpolet of the corresponding order. The conditions (114,115) then factor into d 
separate one-dimensional conditions such that if rnj (q) is taken as the two-scale symbol for an 
/tli-order interpolet in one dimension, mo(q) will be the two-scale symbol for an /^ n -order interpolet 
in d dimensions. Because the two-scale coefficients are the Fourier series coefficients of mo(q), the 
two-scale coefficients for the ci-dimensional interpolet are the outer product of those for the gen- 
erating one-dimensional interpolet, c n = nf=i c n?- Finally, from the two-scale relation (|113|), the 



scaling function corresponding to mo(q) is just the outer product of the one-dimensional generating 
interpolet, T{x) = Yii=x ( x i) ■ As an example, Table |] gives the two-scale coefficients for the 
first-order interpolet in d = 3 dimensions. 

Although the outer-product prescription may always be used to generate higher-dimensional 
interpolets of very simple form, the resulting functions are generally not as compact as possible. To 
illustrate the fact that working in a full three-dimensional framework can generate more compact 
functions, we conclude this section by constructing a third-order interpolet in d = 3 dimensions 
which is not an outer product of three one-dimensional functions. 

For the third-order interpolet in three dimensions, Eqs. (114,115) represent a total of eighty 
linear conditions on the two-scale coefficients. Considering functions with full cubic symmetry 
reduces this to just seven conditions, as we now show. The reflection symmetries about the coordi- 
nate planes, c nijn2in3 = c± nij ± n2i ± n3 , ensure that the first- and third-order conditions, as well as the 
conditions for the off-diagonal second-order multinomials X1X2, X1X3 and X2X3, from (107) are all 
satisfied, leaving only the zeroth-order condition and the three diagonal second-order conditions, 
x\ , x\ and x\. This set of four conditions has implications for mo(q) at q = and the seven points 
irr/i. By the permutation symmetries, among these eight total sets of implications, only those for r\ 
at (000), (001), (011), (111) need be enforced explicitly. The cardinality condition ( |105[) , moreover, 
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implies that one of these sets is redundant, so that we need not impose the conditions at (Oil). 
The zeroth order conditions at the three remaining points appear as the first three conditions in 
(|128j ). For the diagonal second-order conditions, symmetry renders the three conditions at each of 
the two points (000), (111) equivalent. One condition from each of these points appears in (128). 
At the final point, (001), the conditions for x\ and x\ are also equivalent, and so only one of these 



appears in (128). Finally, the x\ condition at (001) makes up the last of the seven conditions in 



A 



c {000} 
c {001} 
c {011} 
C{111} 
c {003} 
c {013} 
V C {H3} J 
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Table [| gives the resulting coefficients along with the coefficients for the outer product of three 
one-dimensional, I = 3 interpolets. Figure |l7] illustrates the more compact nature of the result of 
the three-dimensional construction by comparing the supports of the two functions in the X3 = 
plane. The three-dimensional volume of the support of the new function is smaller than that of the 



product form by a factor of two. Finally, Figure 15 shows the appearance of the new function in 
the X3 = coordinate plane as computed from the recursion (126). 



VI Multilevel Methods 



As Section |V-B.2 discusses, the direct approach of applying the linear operators T, J , Z', j\ 



O, L by multiplying by the corresponding matrices is relatively costly, requiring thousands of 
operations per coefficient in the restricted multiresolution analyses typical in the calculation of 
electronic structure. For unrestricted multiresolution analyses, more methods more efficient than 
direct matrix multiplication exist which exploit the structure imposed by the two-scale relation 
to expend only O(l) floating point operations per expansion coefficient. Many of these methods, 
however, process information on the unrestricted grid and, thus, in a restricted multiresolution 
analysis of n r functions will expend 0(n/n r ) floating point operations per expansion coefficient, 
where n is the number of points in the unrestricted multiresolution analysis. For all-electron 
calculations, n/n r > 10 4 , and these approaches are also extremely inefficient. 

The recent resolution of this problem has been to find new methods which require only 0(1) 
operations per coefficient but which give results which are unchanged when ignoring the n — n r co- 
efficients associated with functions removed in the restriction of the multiresolution analysis. These 
methods then attain the goal of expending only 0{n r ) operations in a restricted multiresolution 
analysis of n r basis functions. We refer to such methods as restrictable. 
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This section reviews the approaches used in the calculation of electronic structure, which in- 



cluded both unrestrictable flWC96| , |TW97[| and restrictable [|BCR9l| , |LAE98| methods. For clarity of 
presentation, we consider in this section only operations in unrestricted multiresolution analyses, 



and dedicate Section VII specifically to the effects of restriction. 



VI- A Transforms: J, Jt, J, jt 

VI- A. 1 General and orthonormal bases 



VI-A.l(a) Forward transform The forward transform of Sec. III-C.l converts the multiscale 



expansion coefficients Fjf-.M into the values of the function on a set of points in real-space. In 
a multiresolution analysis, the product In-.mFn-.m gives the coefficients Fjy of the single-scale 
expansion of the function in Vjv- The values of the function on the points p of the real-space grid 
G are then 



(/)» 



2 N {p-a))(l N:M F N .. M ) 



or in matrix language, 



where 3> 



pa 



(2 N (p — a)). The forward transform operator Z of fli~0|) is thus 



$1. 



N:M- 



(129) 



In general, the invariance of the single-scale basis for Vjy under the translations of the lat- 
tice Cm means that multiplication by the matrix $ takes the form of a convolution. Multiplica- 
tion by <J> may thus be carried out using fast Fourier transform techniques with approximately 
10nlog 2 w floating-point operations |Dau92, WC96|. Alternately, when the <j)[x) are products of 
one-dimensional functions, $ may be factored into one-dimensional convolutions along each of the 
coordinate directions, which expend a total of d(Vol supp 0(x)) 1 / d n operations where d is the di- 
mension of space and Vol supp <f>(x) is the volume of the region of space over which the function 
4>(x) is non-zero. For multiresolution analyses with cardinal scaling functions (both semicardinal 
bases [|ACLT95 , |Yes97 ] and lifted bases |GI98(| ), multiplication by <3? is far simpler: cardinality of the 
scaling functions implies $ = /, where I is the identity. 



Figure [I9j illustrates the flow of information corresponding to each factor Zp+^p in the defini- 
tion ( |sT| ) of Fn-.m- The circles of the upper and lower rows represent coefficients of the multiscale 
expansions Fj^-.p and Fw.p+i, respectively. On the upper row, the larger and smaller open circles 
represent the coefficients for the scaling functions of scale P and detail functions of scale P + 1, 
respectively, and the open circles on the lower row represent coefficients for the scaling functions of 
scale P + 1. The filled circles on both rows represent coefficients from finer scales, which are unaf- 
fected by Tp + i 5 p. The arrows represent the individual multiply-add operations in the application 
of Tp + i 5 p to compute -F/viP+i = ^p+i,pFn-.p- Each arrow multiplies the expansion coefficient at its 
base by a constant factor and accumulates the result onto coefficient at its head. The factors which 
the arrows carry are the two-scale coefficients c n and d n , or unity in the case of the dotted arrows. 
Note that the coefficients Fn-.p+i are taken to be zero before the accumulation begins. Also, for 
clarity, the figure shows only the connections associated with two-scale coefficients in the range 
|n| < 1. Figure 20 illustrates, for three levels, the cascading data flow of ( |30|) when the Tp + i p 
string together to make the final multiscale operator In-.m- 

The translational symmetry of the c n and d n coefficients, reflecting the invariance of the two- 
scale basis for Vp © VFp+i under the translations of Cp, implies that the data flow corresponding to 
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Figure 19: Information flow of the operation F^-.p+i 
ysis). 



1p+\ pFn : p (general multiresolution anal- 



2p+i p takes the form of a set of convolutions, one for each scaling and detail function of the two- 
scale basis. As with these convolutions may be carried out using Fourier transform techniques, 
or, in the case of product basis functions, factored into sets of one dimensional convolutions along 



the coordinate directions and treated all at once[Dau92, Ch. 10] 



VI- A. 1(b) Conjugate forward transform The conjugate forward transform is 

Tt - rt $t 



(130) 



The implementation of this expression parallels that of ( |129| ). In the general case, is a con- 
volution carried out as described above, and for the case of cardinal scaling functions, 3>t = /. 
The component factors F P+1 P of F^ N . M are again invariant under the translations of Cp and may 
be either implemented with either Fourier techniques or, for product functions, factored into one- 
dimensional convolutions. 

The information flow for multiplication by mL.m ^ s a S a i n a cascade of two-scale connections 



1 



t 



but now in reverse order. Note that, when written in terms of components, the matrix 



multiplication a = Mb is a\ = J2j mijbj, and multiplication by the Hermitian conjugate, a = Aft 6, 
is aj = J2i m ijbi- This means that multiplication by the Hermitian conjugate simply reverses 



the direction of and conjugates the factors associated with the arrows. Figure |21| illustrates this 
information flow for one stage I p+1 P of the cascade for F^ N . M - 



VI- A. 1(c) Inverse transform For general unrestricted multiresolution analysis and restricted 
semicardinal multiresolution analyses, the forward transform I is invertible, and we thus consider 
the case J = J -1 . Then, the inverse transform is 



J 



(131) 



For orthonormal wavelets, the inverse convolution must be computed using Fourier techniques 
which require access to the unrestricted grid. The use of cardinal scaling functions removes this 
complication, for then simply = /. 

In multiresolution analysis, the condition ( ffo|) ensures that the two-scale components F-p\\ p 
of the cascade for F~^ l M always exist as linear operators. The invariance of Tp + i 5 p under the 
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Figure 20: Information flow of the multiscale operation Fn : n = In-.mFn-.m (general multiresolution 
analysis) . 
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Figure 21: Information flow of the operation -F/v : p = 2p + i pi*jV:P+i (general multiresolution anal- 
ysis). 
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Figure 22: Information flow of the operation Fjy.p 
ysis). 



2p + i pFN:P+l 



'general multiresolution anal- 



translations of Cp ensures that the information flow for Ip]_i P always takes the form of Figure 22 



for some set of dual coefficients c n and d n , defined by convention with the complex conjugation as 
labelled in the figure. For orthonormal wavelets 1p\i P = ^p + i p so that, comparing Figures 21 
and 22, the two-scale coefficients and their duals are identical. The dual coefficients also have a 



special structure in semicardinal bases, which we shall explore in Sec. VI- A. 2 



VI-A.l(d) 

gate of flml) 



Conjugate inverse transform The remaining transform is the Hermitian conju- 



(132) 



where for any invertible matrix, = (A^ 1 )^ = (A^) -1 . In general, the operation may be 
implemented either by inverting the convolution associated with & with Fourier techniques. Again, 
cardinal scaling functions are much simpler and give just & = I. The operations of the cascade 



for 1 N . M are obtained by complex conjugation and reversal of the flow of Figure E2. As may be 



seen from Figure |23j, the convention for the definition of the dual coefficients has been established 
to ensure that the conjugate inverse transform is just the forward transform associated with the 
dual coefficients. In the orthonormal case, where the coefficients and their duals are the same, this 
represents the expected result that the conjugate inverse and forward transforms are the same. The 
next section discusses the form for the 1p\i p in semicardinal bases. 



VI-A.2 Semicardinal bases 

As we saw in the preceding section, orthonormal bases have the advantage that the dual coeffi- 
cients appearing in the inverse transforms are identical to the two-scale coefficients of the forward 
transform. The disadvantage is that access to the unrestricted grid of finest resolution is needed 
to implement the inverse convolutions <I>~ 1 and <I>~^. Because these later operations are trivial for 
bases with cardinal scaling functions, semicardinal bases will have a great advantage provided the 
dual coefficients needed for the inverse transforms are simple in form, as we now verify. 

The particularly simple form for the dual coefficients for semicardinal multiresolution analyses 
arises from the structure of the Ip + i t p. To elucidate this structure, we analyze Zp+^p into its 
effects on different scales, Xp+^p = (Vc P+1 + Vb f+1 )Zp+1,p('Pc f+1 + T > b p+1 ), where -Bp+i is the 
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Figure 23: Information flow of the operation F^p+i = Zp^ x pFn-.p (general multiresolution anal- 
ysis). 



set of points on scales beyond P + 1 as defined in (p6[). Because Zp + \ t p does not affect functions 
beyond scale P+ 1, we may use (|86|) to simplify this to Zp+^p = Vb p+1 +T > c p+1 ^p+i,pT > c p+1 ■ Next, 
according to (|92|), the fact that the two-scale basis for Vp + Wp+i is semicardinal on the hierarchy 
C P C Cp+i gives V Cp+1 Ip+i,p'Pcp = Vc P + 'P Dp+1 1p+i,p'Pc P and Vc P Ip+i,pV Dp+1 = V Dp+ ^- 
Combining these results, we have 



Zp+i,p = V Bp+1 + {(Vc P +V Dp+1 l P+ i,pVc P )+V Dp+1 ) 
= I + Vd p+1 Zp+i,p'Pcp ■ 



(133) 



As Figure [24|a illustrates, this means that, apart from the simple copy operation /, the only effect 
of Zp + i ; p in the semicardinal case is to add onto the detail points -Dp+i information coming from 
the coarse points Cp. 

The full forward transform, given this structure for Zp + i 5 p and the fact that $ = J, therefore is 



1= J] {i + r Dp+1 i P+1 ,pV Cp ), 



P=N-1 

which corresponds to the following recursive algorithm for computing / = IF, 

Fn-.m = F; 

Fn-.P+I = Fn : p +Vd p+1 ^P+1,pT , C p Fn:P; 
f = Fn-.N- 



(134) 



(135) 



As we show in Sec. [VII , this algorithm is restrictable and gives the same results even when com- 
pletely ignoring points and coefficients eliminated in a restricted multiresolution analysis. Thus, 
with this approach, the forward transform is extremely efficient and may be computed with only 
Vol supp Z(x), or d(Vol supp l(x)) 1 / d in the case of a product interpolet, operations per coefficient. 
As described above, the conjugate to the forward transform (|130|) must simply reverse the 



direction of the information flow relative to Figure 24 a. In product form, it is 



7V-1 



x] = n (i+-p CP ip +1 , P r DP+1 ), 



(136) 



P=M 
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Figure 24: Information flow in a semicardinal basis: (a) forward transform, (b) inverse transform. 
(For respective conjugate transforms, conjugate and reverse the direction of the arrows.) 
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corresponding to the following algorithm for computing F =X*f, 

Fn-.n = /; 

FN:P = FN:P+l+T > Cp1p+\,pPDp +1 FN:P+l', 
F = Fn-.M- 



(137) 



Despite the fact that information appears to flow in from missing points when reversing the arrows 



in Figure 24a, we show in the next section that this algorithm is also restrictable and thus also 
efficient in restricted multiresolution analyses. 

Determining the dual coefficients needed to invert each stage Ip+i p amounts to finding a 



prescription for undoing the effects of the forward transform. As Figure 24a shows, semicardinality 



ensures that the forward transform has no effect on the data sitting on the points of Cp, so that in 
the inverse transform too, these data may be left unmodified. The information flowing diagonally 
from the points of Cp in the forward transform, however, does corrupt the values on the points 
of Dp. Because the information which flowed along these links sits on the points of Cp at the 
start of the inverse transform, to determine the values originally sitting on Dp, one need only 



carry information along the same diagonal links but with coefficients the opposite sign. Figure 24 b 
illustrates the resulting data flow for the inverse transform. Note that Figures p4[a-b are constructed 
so that the entire diagram may be viewed from top to bottom as the process which first performs 
the forward transform and then inverts it. Note that, to accomplish this, the sequence of spaces in 



Figure 24 b is reversed relative to that of the preceding diagrams. 

l semicardi 

c n . Expressed in our matrix language, Figure 



Comparing Figure 24b with Figure 22 we see that in semicardinal multiresolution analyses, the 
dual coefficients are c* = d r , 
24 b represents 



6 n o and dt 



2d r , 



-i 



-p+i,p 



I - 7 y D P+1 Tp+i,p'PCf 



(138) 



which may be verified algebraically by multiplication with ( |133|) and use of the identity VcpVDp+ x 
0. The full inverse transform is thus 

N-l 



j=t- x = n (i-v Dp+1 ip +1 ,pv C p), 

P=M 

which may be cast into a recursion to compute F = J f: 

Fn-.n = /; 

Fn-.P = Fn-.P+i — Vd p+1 Ip+1,pPc p Fn:P+1] 

F = F N:M . 



(139) 



(140) 



Section |VII verifies that this inverse transform is also restrictable. 
Finally, the conjugate inverse transform is 

M 

^= II {l-Fc P lp +1 ,pVn P+1 ), 

P=N-1 

which may be computed with the restrictable algorithm, 

Fn-.m = F; 
Fn-p+i = Fn : p — Vc P Zp + i t pT > D P+1 FN:p; 
f = F N:N , 



(141) 



(142) 



which simply reverses the direction of information flow from that of Figure p4|b. 
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Figure 25: Information flow for a lifted semicardinal basis: (a) forward transform, (b) inverse 
transform. (For respective conjugate transforms, conjugate and reverse the direction of the arrows.) 
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VI-A.3 Lifted bases 



We close this section with a brief consideration of bases constructed by lifting semicardinal bases 
through the scheme of [3we96] as in GI98|. In the lifting scheme, each detail function ip a (x) 



consists of an original un-lifted detail function ipa\x) plus an arbitrary linear combination of 
scaling functions from the next coarser scale, 



ip a (x - r/c/2) = tp^H x - rj a /2) + s a:n <p(x 



n) 



where the r\ a are the decoration points on which the detail functions are centered. (See Eq. (ff"4|).) 
Although a non-zero choice for the coefficients s Qjn destroys semicardinality, these coefficients may 
be chosen to optimize the basis for other purposes. In [GI98], for instance, this combination is 



taken to produce detail functions with zero integral and zero multipole moments up to some order 
for use in solving Poisson's equation. 

With detail functions of this form, each stage of the cascade for the full forward transform may 
then proceed in two phases: (1) the lifted multilevel description i*V:P of the function is converted 
to the un-lifted representation F^} p by transferring the weights s aj7l from the detail functions of 
scale P + 1 to the scaling functions of scale P; (2) the un-lifted forward transform is carried out 



on F^} p to produce FjV:P+i- Figure 25 a shows this two-stage process for the component operation 
2p+i p in the forward transform of a lifted a semicardinal basis. The conjugate forward transform 
just reverses this information flow and conjugates the factors carried by the arrows. 

One may apply precisely the same logic as used in generating Figure p4| b to invert each of these 
two stages in sequence to produce the lifted version of 1p\i P - Figure p4| b shows the resulting 
information flow. The conjugate to this inverse reverses the information flow and conjugates the 
factors carried by the arrows in the figure. 

Finally, we note that to compute the full transforms, one cascades the above component opera- 
tions according to (129- 132) ) and uses the fact that, because the scaling functions remain the same, 



<3? = / is still true. 

VI-B Operators: O, C 

When working with unrestricted multiresolution analyses, the operation of multiplying by matrices 
of the form 

M a p = j d d xb* a {x)Mb p {x), (143) 

where M is some integral-differential operator, is most easily performed in the finest single-scale 
representation (V/v) where the matrix elements are invariant under translations of the lattice Cat 
and thus represent a simple convolution. Substituting the construction formula 

peC N 

into (|143|) gives the overlaps in the multiresolution representation as 



M = X\, M MX N ,M, (144) 

where M is the single-scale overlap matrix M on the finest scale defined in ( |116| ) of Sec. [V-C.2] . 
This immediately gives an 0(n) method for applying A4, where n is the number of points on the 
finest grid Cn'- transform to the single-scale representation on level N, apply the convolution M, 
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and transform back to the multiresolution representation with X*. The main disadvantages of this 



approach, which is used in [WC96, TW97], are that (a) it is not restrictable and (b) it requires 
storage and processing of data on the unrestricted grid Cjv- 

The non-standard multiply approach of Beylkin, Coifman and Rokhlin[BCR91] provides one way 



of applying operators in multiresolution bases with linear computational effort while circumventing 
the need to work directly with the Vn representation. A more general approach [ L AE9§| ] , which 



includes the non-standard multiply as one special case, is to follow the strategy which underlies 
Greengard and Rokhlin's fast multipole method [pR8S |. The idea is to eliminate the redundancies 



inherent in the matrix elements of Ai by lumping information together into coarse and fine con- 
tributions. Only, rather than using multipole moments to summarize fine-scale information in an 
approximate way, one may exploit the two-scale relation to gather together the fine-scale overlaps 
exactly. And, rather than using Taylor series to summarize smooth information approximately, 
one may use the two-scale relation to represent smooth-scale information exactly in terms of coarse 
scaling functions. The remaining information on proximate scales then may be handled efficiently 
by direct, and thus also exact, multiplication. 

To define precisely the separation into "fine," "proximate," and "smooth" contributions, we 
first decompose the operator Ai into contributions affecting the results on the different scales 
M < Q < N of the multiresolution analysis, 

M = (V Cm + ... + V Dq +... + V Dn )M. (145) 

For a given term specified by its scale Q, there are contributions which originate from scales R 
which are either "proximate" (\R — Q\ < £), "finer" (R > Q + £) or smoother (R < Q — £), where i 
is a measure of scale separation. Following this classification, for a given Q, we partition the set of 
all points in the multiresolution analysis as 

C n = SqUPqUFq. (146) 

Although the interpretation of Sq, Pq and Fq is clear, the precise nature of these sets of points 
involves a variety of minor details as Q approaches the coarse and fine limits of the multiresolution 
analysis. In particular, when Q — I < M, all of the smoother grids (including Cm) are proximate 
and Sq is empty. On the other hand, when Q + 1 > N , all of the finer contributions are proximate 
and Fq is empty. An additional notational complication is the fact that the points on scale Q = M 
are the lattice Cq whereas the points on the finer scales Q > M make up the detail crystals Dq. 
To ease the burden of constantly making this distinction, we define 



f C Q if Q = M 
~ D Q if Q>M 



as all points on the "level" Q. With these considerations, the precise definitions for Sq, Pq and 
Fq are 



Sq = 
Fq = 



if Q - £ < M 
C Q . e if Q-£>M 

if Q + £ > N 

Bq+i-! \iQ + £<N 



min(N,Q+£-l) 
Pq = C N -Sq-Fq= (J Lr, 

R= max(M,Q-£)+i 
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where is the empty set. With these definitions, the full decomposition of A4 is 



A' 



M = £ V Lq M(V Sq + V Pq + V Fq ) (147) 



Our discussion for the implementation of these various terms will be limited to the semicardinal 
case, which has been shown to be restrictable|LAE98], and so we return to the specialized notation 



T(x) for the scaling functions, in place of the generic <fi(x). 

To evaluate the proximate contributions we first define the inter-scale matrices of overlaps 
connecting the scaling functions of Vp with the scaling functions of Vr, 

f fd d xl*(2 p (x-p))Ml(2 R (x-r)) if peCp and r£ C R , . 

(Mp >^" otherwise ' (148) 



which, following the same argument which led to (144), take the form 

Mp,r = Vc p 1 ] n ., p M1 n ..rVc r ■ (149) 

Note that the Mp t p are invariant under the translations of the lattice C m i n (p^) an d thus may be 
implemented as convolutions. 

The proximate contributions all involve terms of the form VlqMVl r = VlqZ* MTVl r - To 
simplify this expression, we use the identity 

Tn-.m'PLr = In:rPl r , (150) 



which may be proven using (31,133). Alternately, Eq. ( |150| ) may been seen immediately as the 
algebraic statement that in a semicardinal multiresolution analysis the basis functions associated 
with the points p of scale R, whose values are the columns picked out of In-.m by Vl r on the left of 
( |15C| ), are just the scaling functions from the single-scale representation of scale R associated with 
the same points p. With ( |150() , the proximate scale terms are just 



Vl q MV Lr = V Lq iI m M1 N:M V Lr (151) 
= V Lq V Cq 1 ] n .. q M1 n ..rVcrV Lr 
= V Lq M q ,rV Lr , 

where on both sides of M in the second line we have used ( |15C| ) and the fact that always Lr C Cr 
so that V L r = VcrVlr- 

Unlike the stages of the transforms, the inter-scale matrices here used to apply operators do 
not convert one multiscale representation into another but rather link together different single-scale 
representations. We thus introduce a second useful diagrammatic representation for the coefficients 
of a multiresolution analysis. The horizontal rows of Figure ^ represent the scaling functions for 
each single-scale representation Vq. The basis for a semicardinal multiresolution analysis consists 
of a subset of these functions and includes for each point only the scaling function of the coarsest 
scale containing that point. These scaling functions are indicated in the figure as filled circles. In 
practical implementations it is convenient to maintain data structures which include coefficients 
for the redundant scaling functions from finer scales, which are represented in the figure as open 
circles. We refer to this as the redundant representation. This redundant representation simplifies 
significantly the implementation of the Mq^r as convolutions and represents only a small increase 
in storage by a factor of 2 d /(2 d — 1), which is just 8/7 in d = 3 dimensions. (See [ LAE98| .) 
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Figure 26: Information flow for proximate contributions to Ai: inter-scale convolutions Mrq 
(arrows). 



The arrows in Figure 26 represent the calculation of the proximate contributions according to 



(147 ,151). The calculation begins with the coefficients of the vector to which the operator is applied 
residing on the left of the diagram and with all coefficients residing on the right being initialized 
to zero. Then, each arrow pointing from a level R on the left to a level Q on the right accumulates 
onto level Q the result of the convolution Mq^r applied to the information present on level R. Note 



that the factors Vl r in (151) imply that the coefficients associated the redundant functions on the 
left of the figure should be set to zero before performing these convolutions. After the accumulation 
of the convolutions, the proximate contribution to Ai appears on the filled circles on the right side 
of the diagram. As with all diagrams in this section, the factors Vl q in ( |151|) indicate that the final 
values computed on the empty circles are extraneous and simply may be ignored, or not computed 
at all. 

Next, we turn to the smooth contributions, which need only be computed when Q — i > M. 
The basic strategy for these terms is to expand the smooth contributions exactly in terms of scaling 
functions on scale Q — i. To do this, note that with the expression ( |I44| ) for A\ and the definition 
Sq = Cn_£, Vd M.Vs ends in a product of the form In-.mVsq = In-.mVcq which becomes 



J Q-t ' 



-JV:t 



*XPc Q . 



-N:M ) 



(152) 



after some manipulation of ( |134| ), The result (|152j ) also may be seen immediately as the algebraic 
statement that in a semicardinal multiresolution analysis the multiscale basis functions from scale 
Q — £ or coarser, whose values appear in In-.m'Pcq-^ niay be expanded exactly in the single-scale 
representation on level Q — t with coefficients equal to the values of these basis functions at the 
corresponding points of Cq-h. Using (|I50|) on the left, (152) on the right, and the result ( |I49|) , the 
smooth contributions become 



Vl q MVs q = V Lq 1 ] N:Q M1 n .. q _{P Cq _ 1 1 

= Vlq^Q,Q-(Xn:M- 



N:M 



(153) 



Figure illustrates the information flow for the smooth contributions ( |I53j ). First, the forward 



transform 2n-.m is carried out using the cascade algorithm (135), as represented by the curved 
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Figure 27: Information flow for smooth contributions to ftA: cascade algorithm for In-.m (curved 
arrows), access to data on multiple scales (vertical arrows), inter-scale convolutions Mr^q (diagonal 
arrows) . 



arrows. (In practice, the cascade need not be carried out completely as the results on the finest 
few levels will not be used.) Next, the vertical arrows indicate that one must ensure that the 
resulting coefficients are replicated on all finer levels, for they will be accessed also from these levels 
by the Mq^q-i. The final step in computing the smooth contribution is to perform the inter-scale 
convolutions Mq q_£ represented by the diagonal arrows as in Figure |2^. 

The strategy for the finer contributions, present whenever Q — i < N, is to gather overlap 
information from finer scales in the form of the matrix elements, 

J"q=V Cq MV Fq . (154) 



Note that the finer scale contributions in (147) for any scale Q may be extracted from these objects 
as simply VlqFq. To gather the effects of finer scales on these matrix elements, we use the 
decomposition Fq = D Q+i U F Q+ i and make the replacement 1n-.qVc q = Zn:Q+iPc q+1 Zq+i,q'Pcq 
to convert (|154j) into the recursion 



FQ = Vc q MV Dq+1 +V Cq MV Fq+1 (155) 
= m q,q+iPd q+1 + Vc q Zq +1iQ Fq+i , 
which may be solved for !Fq with T^-i = 7 ) Cq-M-'Pd n as the starting point of the recursion. 



Figure |2q illustrates the information flow for this process. The Mq^q + i convolutions may be first 
computed and then "folded- in" as the recursion (|155| ) proceeds. These convolutions are indicated 
by the diagonal arrows. As with the proximate contributions, the projections Vdq +1 imply that the 
redundant points associated with the open circles on the left carry zero value for these convolutions. 
Once the convolutions are precomputed in this manner, the finer scale contributions are computed 
from the recursive accumulation of the action of q from each scale onto the next coarser scale 
in sequence, from the bottom to the top of the figure, as represented by the curved arrows. In 
practice, the operations Fq +1 q on the finest few scales work with data with zero value and so need 
not be computed. 
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Figure 28: Information flow for fine contributions to Ai: inter-scale convolutions M^q (diagonal 
arrows), 2p +1 p stages proceeding in sequence from finest to coarsest (curved arrows). 



The sum of the contributions illustrated in Figures |2 6^28] gives the final result of the application 
of the operator M. The complete procedure for computing H = MG for any vector G for the 
non-trivial cases N — M > £ which involve more than just proximate contributions is thus 



N-e 



M 



N—£,N 'Dff 



Vd n G, g = lG; 



(156) 



Tq = M Q> Q + eVD Q+e G + V Cq 1 ] q+1 ^q+i (for Q < N -£); 



N 

q=m+£ 

N / mm(N,Q+l-i) 

+ E E v Lq m q>r v Lr 

Q=M \r= max (M,Q-£+l) 
N-e 

+ E Tl^q. 

Q=M 



The non-standard multiply approach of Beylkin, Coifman and Rokhlin[BCR91] corresponds to 
the special case i = 1 case of this approach. For the restrictions usually encountered in atomic 
calculations, the case £ = 2 has been shown to improve the speed of calculations by at least a factor 



of four[LAE98, Lip9S]. 



VII Impact of Restriction 

We now consider how the prescriptions ( 135| , 137 ,14C , 142 ,156) for operations in semicardinal mul- 
tiresolution analyses are affected when the analysis is restricted. We show that under relatively 
mild conditions on how quickly the restriction changes resolutions, it is possible to simply ignore 
the missing coefficients and yet obtain exactly the results of an unrestricted multiresolution analysis 
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and thus achieve the ultimate goal of expending O(l) operations per expansion coefficient while, 
remarkably, introducing no additional approximations beyond the choice of basis set. These issues 



were first explored in [LAE98], which prompted the study and development of the concept of "S- 
trees" to appear in a coming work[ |CD97 |. We now present the discussion of [LAE98] using the 
explicit matrix language introduced in Sec. IV-F . 

We begin by introducing the notation [*] as the restriction of the scope of the object "*" to the 
functions and grid points surviving the restriction of a multiresolution analysis. The restriction [S] 
of a set of points S, for example, is defined as the set of those points of S associated with functions 
which survive the restriction. The set of all surviving points is the restriction of the finest grid in 
the multiresolution analysis, [Cat]. Thus, in general, for any set of points S, 



[S] = [c N ] n s, 



(157) 



or in terms of projections, 

T\s\ = V [C N ]Ps 

The restriction of any linear operator O is defined so that its matrix elements ignore coefficients 
that have been restricted from the multiresolution analysis. Algebraically, this is achieved by 
projection, 

[0] = V [ c n] OV [ c n] . (158) 



Consistent with this definition, the matrices O and L used in Sec. III-C.l involve only overlaps 
between functions present in the basis. 

This simple prescription for the restriction of linear operators, however, is not appropriate in 
general for the transforms, which relate expansion coefficients on one side with function values at 
points in real space on the other. Even a smooth function will have significant values in real space on 
the grid points associated with extremely fine functions. Thus, great care is needed to ensure that 
the lack of access to this information does not adversely affect the results of expressions involving 
the transforms. For semicardinal bases, as we shall see, very mild conditions on the restriction of 
the multiresolution analysis can ensure that the unknown expansion coefficients will have absolutely 
no impact on the function values computed by the forward transform on the surviving points, so 
that 



(159) 



and, most critically, that the missing real-space information has absolutely no impact on the ex- 
pansion coefficients computed by the inverse transform for the surviving functions, 



[J] = V [Cn] JVi 



[CV 



(160) 



Taken together, conditions ( |159| , 160| ) have the remarkable consequence, to which we first alluded 
in Sec. Q-B , that, even in a restricted basis, we may compute the expansion coefficients for any 
non- linear combination of fields, for example the exchange-correlation energy density e xc (r), and 
obtain results identical to what would be obtained by working with an unrestricted multiresolution 
analysis of arbitrarily high resolution. 



VII-A Transforms 

For the transforms, we begin by considering the implications of condition Q16C| ) for the allowable 
restrictions of semicardinal multiresolution analyses. The inverse transform in semicardinal bases 
(|13S| ) has a particularly simple form which is useful for understanding the impact of restriction. 
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Because Vcp'Pdq+i = for all P < Q, the fact that the product (139) orders terms from coarsest 
to finest from left to right, implies that all but the zero and first order terms in T > d p+ iIp+i,p'Pcp 
vanish from the product, leaving only the sum 

N-l 

J = I — /, T > D P+l 1p+\,p'Pcp- 



P=M 



Inserting this into ( 160 ) and using ( |157| ) gives 



N-l 

E n 

P=M 



N-l 

V [Cn] ~ E ' P [Dp +1 ]Zp+1,p'PCf 
P=M 



Because the [Dp+J are disjoint, the sums must be equal term by term in V[d [ 
condition 



T , [Dp +1 }3p+i ) p{'Pcp ~ V 



V[D P+1 fp+l,pV C p-[C P } = 0. 



leading to the 
(161) 



Given that all the points in [Cp] must appear in [Cp+i], the simplest geometric interpretation of 
( |161| ) is that no single-scale basis function associated with a point dropped from scale P may have 
a non-zero two-scale expansion coefficient for any scaling function maintained in the basis on level 
P + 1. This is referred to as the good- grid condition. 

Figure 29 illustrates the effects of the good-grid condition on the progress of a calculation. The 
figure follows the same conventions as the figures of Sec. VI- A , but now the multiresolution analysis 
is restricted to provide the resolution of Vm in the left half of the figure, to make a brief transition 
through Vm+i, arid finally to attain the resolution of Vm+2 = Vn on the right. Although the figure 
shows the information flow which would occur in the full multiresolution analysis, it shows only 
those coefficients which survive the restriction. Geometrically, the condition (161) states that no 
arrow in the forward transform (upper half of the figure) which originates from a missing point 
may terminate on a surviving point. This ensures that the surviving coefficients computed at 
each stage of the restricted transform are unaffected by the missing coefficients. For semicardinal 
multiresolution analyses, the pattern of the information flow for the inverse transform (lower half of 
the figure) is the same as for the forward transform. Thus, the inverse too is unaffected by ignoring 
coefficients on the restricted points. Finally, the figure illustrates the observation made above that 
so long as the good-grid condition is maintained, the process of computing expansion coefficients 
for the results of a local non-linear interaction (represented by the vertical solid connections in the 
center of the figure) such as e xc (J2i fl^i^)] 2 ) remains unaffected by coefficients restricted from the 
multiresolution analysis. 

Figure illustrates the constraints of the good grid condition on allowable restrictions for 
functions of support ±2. The presence of the one point on scale M + 2 in the figure requires the 
survival, in the restriction of the multiresolution analysis Vm © Wm+i © Wm+2, of all of the other 
points in the figure. As the figure illustrates, the physical requirement of the good-grid condition is 
just that the restriction pass through finite regions representing each level of resolution in sequence 
and that it not not skip discontinuously from level to level. Because the expansion coefficients 
of physical functions exhibit this behavior naturally (Figure ^), the good-grid condition generally 
represents little or no additional burden in the restriction of the basis. 

To establish algebraically that on good grids the procedures ( p35|Jl4C|Jl3^JT42[ ) give correct 
results when working only with data on the surviving points, we first note that on a good grid the 
separate stages of the transforms individually obey the restriction condition, 



p+i,p 



l,p 



1 



-1 

p+l,p 



-1 



P[C N ]lp+i,p, 



(162) 
(163) 
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Figure 29: Information flow in the calculation of expansion coefficients for e xc in a restricted 
semicardinal multiresolution analysis on a good grid: forward transform (upper half), non-linear 
local interaction (vertical connections in center of figure), inverse transform (lower half). 
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Figure 30: Implications of good grid condition on scales (M : M + 2) stemming from a single point 



(solid circle) on scale M+2 for functions of support ±2: requirements of V^p I 
(solid arrows), points already required from lower levels (dashed arrows). 



_ 1 ]Zp+i,p'Pcp-ICp} = 



as follows directly from the defining good-grid condition ( |161| ) and the decompositions ( 133| , 138| ), 
Next we note that working with only data on the surviving points amounts to ignoring at each 
step of a procedure the input from and output onto the restricted points. Algebraically, this 
corresponds to replacing the factors in the products making up the transforms ( 134|jl39 , 136| , |14l| ) 
with their restricted counterparts, [Zp+^p], [Xp^_ ip ], [Xp +ip ], [Ip+i p] , respectively. 



From these considerations, we see that the forward transform 
coefficients gives 



executed with only surviving 



M 

P=N-1 



] - *P[c n \£n,n-iP\c n \ ■ ■ ■ 'P[c n ]Zp+2,p+iP[c n ] i p+i,p'P[c n ] 

= 'P[c n ]Zn,n-i'P[c n ] ■ ■ ■ 'P[c n ] :i: p+2,p+i {'P[c n ]Zp+i,p'P[c n ] ) 
= 'P[c n ]Zn,n-iP[c n ] ■ ■ ■ (P[c n }Zp+2,p+i'Pic n ])Zp+i,p 

= r [CN] i, (164) 

where we have collapsed the product telescopically using (|l62j ). Directly analogous considerations 
for the inverse transform lead to 



N-l 

n 

P=M 



3-P+l,p\ 



(165) 



Thus, we see that applying the algorithms (135,140) to only the surviving coefficients leads to 
results on the surviving points which are identical to what would be obtained with the unrestricted 
transforms I and J . Note that because the left-hand sides of (164, 163 ) are unchanged by post- 
multiplication by V\p N \i these results also confirm directly that ( |15£ , 160[ ) obtain on good grids. 

Finally, we note that the conjugate transforms defined in Sec. fl [-E| for any basis are the Her- 
mitian conjugates of the associated forward and inverse transforms. From ( |159 , 160| ) and (164, 165| ), 
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these are thus simply 

N-l 

= n i4 + i,p] 

P=M 



P=N-1 



Composed entirely of restricted component operations, these expressions are precisely what the 
algorithms ( 137|J142 ) compute when working only with data on surviving points. 



VII-B Operators 



The procedure ( |156j ) for applying physical operators divides the contributions to M into a sum of 
three classes of contribution: smooth, proximate, and fine. The proximate contributions call upon 
Mq p to collect data from functions present in the basis onto coefficients present in the basis. Thus, 
regardless of the form of the restriction, the restriction of these parts of M amounts to ignoring 
coefficients restricted from the basis, which is precisely what is accomplished when carrying out 
( |156| ) on data associated with only the surviving points. Algebraically, this corresponds to the fact 
that restricting ( |151| ) leads to corresponding expressions with all component operations replaced by 
their restricted counterparts. The proximate contributions in ( |156| ) are therefore always restrictable. 

The smooth contributions, on the other hand, do impose conditions on the restriction. The 
contributions (153) connect the single-scale basis associated with Cq_£ through the operator A4 



onto the multiscale basis functions associated with the points of Lq. To ensure that no information 
is lost when computing these contributions, the restriction must be such that all scaling functions 
associated with the lattice of the scale i levels coarser which overlap with, or touch, a function in 
the basis must also appear in the basis, 

-PlL Q ]M Q , Q -iV CQ _ e _ [CQ ^ = 0. (166) 

Algebraically, this condition allows the transition from the first to the second line in the derivation 

[PL Q M Qt Q-ll N .. M ] = r [ L Q ]MQ > Q„iPc Q _ £ lN:M'P[C N } (167) 
= T> [Lq] M Q,Q-(P[C q ^ 1 \^N:M'P [ c n ] 

= [Pl q }[M q>q ^][I n ..m], 
which establishes that the coarse contributions to [M] may be computed correctly while processing 



data on the surviving points only. Figure |31| illustrates the requirements of the this condition for 
functions of support ±2 when using the £ = 2 version of the algorithm. 

Finally, we turn to the finer contributions. The first term in the recursion (|155|) for these 



contributions is of the same generally restrictable form as the proximate contributions and so 



impose no condition on the restriction. However, proper computation of the second term of (155) 
requires a condition on the grid which is slightly stronger than ( |166j ) . Now we must require that all 
scaling functions from the scales at least £ levels coarser which touch a function in the basis must 
also survive the restriction, 

V Cq -[c q ]MV [Fq] = 0, (168) 



which, from the definition (154) of J 7 , is equivalent to 3~qP\c n ] = [3~q]- Combining the Q + 1 



case of this condition with the previously established good-grid condition (|161|), we see that when 
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Figure 31: Requirements of the grid-touching condition (168) on scale Q which stem from the 
presence of a detail function (filled circle) on scale Q + 2: required points (solid circles), optional 
points (dashed circles), support of functions (braces). (Illustrated for the t = 2-level algorithm.) 



computing on only surviving coefficients, the second term in ( |155|) gives 



[r CQ ][ih +1 nW Q , 



Vu 



(169) 



= r [C Q ] I Q+l,Q :F Q+l' P lC N ] 
= I^CqIq+i^Q+i], 

precisely the correct restricted result. 

We refer to condition ( |168j ), which contains (166) as a special case, as the I- level grid-touching 
condition. Combined with the good-grid condition (161), this condition ensures that the procedure 



given in Sec. VI-B gives correct results for integral-differential operators even in restricted mul- 
tiresolution analyses. Note that the £-level grid-touching condition requires that the the situation 
illustrated in Figure [3l] holds not only for the scaling functions on scale Q but also for all coarser 
scales R < Q as well. 

As the scale separation £ increases, the detail function on the lower scales appear relatively 
smaller, and the grid-touching condition becomes weaker. In the limiting case of very large I, the 
support of the detail functions appears like a single point and the touching condition approaches 
the good-grid condition. For low £, the differences in the conditions on the restriction can have 
important consequences. For example, going from I = 1, which corresponds to the non-standard 
multiply of [BCR91], to i = 2 allows for quicker changes in resolution that have been shown to 
accelerate calculations by over a factor of four for the restrictions typically encountered in atomic 
calculations [ |LAE98| ] . 

As a closing note, the alternative restriction conditions 

'Pc Q _ i -[C Q _ i fNMV[c N ] = 
V [C Q \ L Q+x,Q P c Q+1 -\c Q+i ] = 



would also allow for proper computation of the smooth and fine contributions in ( 167 , 169 ), respec- 
tively. These conditions, however, are "finer-looking" in the sense that each function in the basis 
requires the presence of functions on even finer scales. This leads to an infinite chain of conditions 
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that ultimately requires the multiresolution analysis to be of infinite resolution. As a result, the 
Mevel touching condition is the appropriate condition to impose in real computations. 



VIII Concluding Remarks 

Mallat and Meyer's development of multiresolution analysis provides a tool not only of great use in 
functional analysis and signal processing but also of great potential for the study of physical systems 
exhibiting behavior on multiple length scales. In particular, multiresolution analysis provides the 
first practical possibility for a unified, systematic treatment of core and valence behavior in the 
electronic structure of molecular and condensed-matter systems. The first all-electron density- 
functional calculations of atoms and molecules carried out with this approach ([ ACLT95| , Ari95] 



Sec. ||.) have demonstrated that multiresolution analysis provides an extremely efficient means of 
representing the core and valence electrons simulataneously. Latter work underscores the promising 



possibilities of combining multiresolution analysis with pseudopotential theory[ WC96], developing 



dynamic restriction schemes [JT W97 1 , and using lifted bases for some phases of the calculations [GI98]. 



Finally, the recent development of the theory of fast restrictable algorithms for semicardinal bases 
( jLAE9SH , Sec. [vn]) now paves the way for the first application of multiresolution analysis to the 



calculation of the electronic structure of large, complex systems. 

In closing, the author would like to leave the reader with the understanding that many op- 
portunities remain for making significant contributions to the field of the multiresolution analysis 
of electronic structure. The basic groundwork is now in place, but many interesting and impor- 
tant open questions remain. For example, the differing stiffnesses of the core and valence degrees 
of freedom clearly call for some new form of preconditioning. The significant expense of solving 
Poisson's equation at each electronic iteration indicates that methods which search directly for the 
saddle point of the LDA Lagrangian will be more efficient than present approaches. Also, pro- 
cedures are needed to anticipate accurately the changes in the expansion coefficients of the wave 
functions as the nuclei move. Finally, there is the intriguing possibility that the weakness of the 
coupling between the core and valence degrees of freedom will allow the atomic cores to be treated 
independently and in parallel. 
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A Cofactors of the semicardinal two-scale decomposition matrix 

In this appendix, we determine the cofactors of the matrix 

i^. = (-!)*•»&, 
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where, as defined in the text, the r\i range over the set {0, l} d . 

The central result we use to determine these cofactors is that P 2 = 2 d I, which we show by 
explicit computation, 

k 

d 

d 

= TT (-l)fak)««*)e + fa*)«) 
6=1 {(Vh) e =0,l} 

d 

= + (-l)^)e+(%)< 

e=l 
d 



e=l 

2 5ij. 



With P 2 = 2 d I established, we have P 1 = P/2 rf = (cof P) T /det P, where cofP is the matrix 
of the cofactors of the matrix P. To determine detP, we note that detP = (2 ) , so | detP| = 
{2 d ) . Thus, 



cof P = (det P)(P" T ) = (±(2 d f' 1 ) (^P T ) = ±(2 d f' 1 - 1 P T 



,2 d ' 

Because all entries in the top row of P are unity, the elements of the first column of the cofactor 
matrix, needed at the end of Sec. V-B.3, all have a single constant value, ±(2 d ) 2 -1 . 
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